Changeset 5654

Timestamp:

09/02/10 09:47:07 (15 years ago)

Author:

seroussi

Message:

improved macayeal if coupling

File:

• issm/trunk/src/c/objects/Elements/Penta.cpp (modified) (1 diff)

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

@@ -2607 +2607 @@
 void Penta::CreateKMatrixDiagnosticMacAyeal( Mat Kgg){
         
+        /* local declarations */
+        int             i,j;
+
+        /* node data: */
+        const int    numgrids3d=6;
+        const int    numgrids2d=3;
+        const int    numdof2d=2*numgrids2d;
+        double       xyz_list[numgrids3d][3];
+        int*         doflist=NULL;
+
+        /* 3d gaussian points: */
+        int     ig;
+        GaussPenta *gauss=NULL;
+        GaussTria  *gauss_tria=NULL;
+
+        /* material data: */
+        double viscosity; //viscosity
+        double oldviscosity; //viscosity
+        double newviscosity; //viscosity
+
+        /* strain rate: */
+        double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+        double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+
+        /* matrices: */
+        double B[3][numdof2d];
+        double Bprime[3][numdof2d];
+        double L[2][numdof2d];
+        double D[3][3]={0.0};            // material matrix, simple scalar matrix.
+        double D_scalar;
+        double DL[2][2]={0.0}; //for basal drag
+        double DL_scalar;
+
+        /* local element matrices: */
+        double Ke_gg[numdof2d][numdof2d]={0.0}; //local element stiffness matrix 
+        double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point.
+        double Jdet;
+
+        /*slope: */
+        double  slope[2]={0.0};
+        double  slope_magnitude;
+
+        /*friction: */
+        double  alpha2_list[3];
+        double  alpha2;
+
+        double MAXSLOPE=.06; // 6 %
+        double MOUNTAINKEXPONENT=10;
+
+        /*parameters: */
+        double viscosity_overshoot;
+
         /*Collapsed formulation: */
         Tria*  tria=NULL;
+        Penta* pentabase=NULL;
 
         /*inputs: */
         bool onwater;
         bool onbed;
-
+        bool shelf;
+        int  approximation;
+        Input* vx_input=NULL;
+        Input* vy_input=NULL;
+        Input* vxold_input=NULL;
+        Input* vyold_input=NULL;
+
+        inputs->GetParameterValue(&approximation,ApproximationEnum);
         inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
 
         /*If on water, skip stiffness: */
-        if(onwater)return;
-
-        /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the 
-          bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 
-          the stiffness matrix. */
-
-        if (onbed==0){
-                /*This element should be collapsed, but this element is not on the bedrock, therefore all its 
-                 * dofs have already been frozen! Do nothing: */
-                return;
-        }
-        else if (onbed==1){
-
-                /*This element should be collapsed into a tria element at its base. Create this tria element, 
-                 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
-
-                /*Depth Averaging B*/
-                this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
-
-                /*Call Tria function*/
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreateKMatrix(Kgg);
-                delete tria->matice; delete tria;
-
-                /*Delete B averaged*/
-                this->matice->inputs->DeleteInput(RheologyBbarEnum);
-
-                return;
+        if(approximation==MacAyealApproximationEnum){
+                if(onwater)return;
+
+                /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the 
+                  bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 
+                  the stiffness matrix. */
+
+                if (onbed==0){
+                        /*This element should be collapsed, but this element is not on the bedrock, therefore all its 
+                         * dofs have already been frozen! Do nothing: */
+                        return;
+                }
+                else if (onbed==1){
+
+                        /*This element should be collapsed into a tria element at its base. Create this tria element, 
+                         *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
+
+                        /*Depth Averaging B*/
+                        this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
+
+                        /*Call Tria function*/
+                        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                        tria->CreateKMatrix(Kgg);
+                        delete tria->matice; delete tria;
+
+                        /*Delete B averaged*/
+                        this->matice->inputs->DeleteInput(RheologyBbarEnum);
+
+                        return;
+                }
+        }
+        else if(approximation==MacAyealPattynApproximationEnum){
+                /*retrieve some parameters: */
+                this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
+
+                /*If on water, skip stiffness: */
+                if(onwater)return;
+
+                /*Find penta on bed as this is a macayeal elements: */
+                pentabase=GetBasalElement();
+                tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+
+                /* Get node coordinates and dof list: */
+                GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids3d);
+                tria->GetDofList(&doflist,MacAyealApproximationEnum);  //Pattyn dof list
+
+                /*Retrieve all inputs we will be needing: */
+                vx_input=inputs->GetInput(VxEnum);
+                vy_input=inputs->GetInput(VyEnum);
+                vxold_input=inputs->GetInput(VxOldEnum);
+                vyold_input=inputs->GetInput(VyOldEnum);
+
+                /* Start  looping on the number of gaussian points: */
+                gauss=new GaussPenta(5,5);
+                gauss_tria=new GaussTria();
+                for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                        gauss->GaussPoint(ig);
+                        gauss->SynchronizeGaussTria(gauss_tria);
+
+                        /*Get strain rate from velocity: */
+                        this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+                        this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
+
+                        /*Get viscosity: */
+                        matice->GetViscosity3d(&viscosity, &epsilon[0]);
+                        matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
+
+                        /*Get B and Bprime matrices: */
+                        tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria);
+                        tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
+
+                        /* Get Jacobian determinant: */
+                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+
+                        /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 
+                          onto this scalar matrix, so that we win some computational time: */
+                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+                        D_scalar=2*newviscosity*gauss->weight*Jdet;
+                        for (i=0;i<3;i++) D[i][i]=D_scalar;
+
+                        /*  Do the triple product tB*D*Bprime: */
+                        TripleMultiply( &B[0][0],3,numdof2d,1,
+                                                &D[0][0],3,3,0,
+                                                &Bprime[0][0],3,numdof2d,0,
+                                                &Ke_gg_gaussian[0][0],0);
+
+                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
+                        for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+                }
+
+                /*Add Ke_gg to global matrix Kgg: */
+                MatSetValues(Kgg,numdof2d,doflist,numdof2d,doflist,(const double*)Ke_gg,ADD_VALUES);
+
+                //Deal with 2d friction at the bedrock interface
+                if((onbed && !shelf)){
+                        /*Build a tria element using the 3 grids of the base of the penta. Then use 
+                         * the tria functionality to build a friction stiffness matrix on these 3
+                         * grids: */
+
+                        tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg);
+                }
+
+                /*Clean up and return*/
+                delete tria->matice;
+                delete tria;
+                delete gauss_tria;
+                delete gauss;
+                xfree((void**)&doflist);
         }
 }