Changeset 5203

Timestamp:

08/12/10 14:32:30 (15 years ago)

Author:

seroussi

Message:

Kmatrix coupling

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (6 diffs)
• Penta.h (modified) (1 diff)
• PentaRef.cpp (modified) (1 diff)
• PentaRef.h (modified) (1 diff)

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

@@ -755 +755 @@
                         return;
                 }
+                else if(approximation==MacAyealPattynApproximationEnum){
+                        ISSMERROR("coupling MacAyeal/Pattyn not implemented yet");
+                        CreateKMatrixDiagnosticMacAyeal( Kgg);
+                        CreateKMatrixDiagnosticPattyn( Kgg);
+                        CreateKMatrixCouplingMacAyealPattyn( Kgg);
+                }
                 else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
         }
@@ -827 +833 @@
                 else if(approximation==HutterApproximationEnum){
                         return;
+                }
+                else if(approximation==MacAyealPattynApproximationEnum){
+                        CreatePVectorDiagnosticMacAyeal( pg);
+                        CreatePVectorDiagnosticPattyn( pg);
                 }
                 else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
@@ -2083 +2093 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattyn{{{1*/
+void Penta::CreateKMatrixCouplingMacAyealPattyn( Mat Kgg){
+
+        /* local declarations */
+        int             i,j;
+
+        /* node data: */
+        const int    numgrids=6; // Pattyn numgrids
+        const int    numdof=2*numgrids;
+        const int    numgrids2=3; //MacAyeal numgrids
+        const int    numdof2=2*numgrids2;
+        double       xyz_list[numgrids][3];
+        int*         doflist=NULL;
+        int*         doflist2=NULL;
+
+        /* 3d gaussian points: */
+        int     num_gauss,ig;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double* fourth_gauss_vert_coord  =  NULL;
+        double* area_gauss_weights           =  NULL;
+        double* vert_gauss_weights           =  NULL;
+        int     ig1,ig2;
+        double  gauss_weight1,gauss_weight2;
+        double  gauss_coord[4];
+        int     order_area_gauss;
+        int     num_vert_gauss;
+        int     num_area_gauss;
+        double  gauss_weight;
+
+        /* 2d gaussian point: */
+        int     num_gauss2d;
+        double* first_gauss_area_coord2d  =  NULL;
+        double* second_gauss_area_coord2d =  NULL;
+        double* third_gauss_area_coord2d  =  NULL;
+        double* gauss_weights2d=NULL;
+        double  gauss_l1l2l3[3];
+
+        /* 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][numdof];
+        double Bprime[3][numdof2];
+        double L[2][numdof];
+        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[numdof][numdof2]={0.0}; //local element stiffness matrix 
+        double Ke_gg_gaussian[numdof][numdof2]; //stiffness matrix evaluated at the gaussian point.
+        double Jdet;
+
+        /*friction: */
+        double  alpha2_list[3];
+        double  alpha2;
+
+        /*parameters: */
+        double viscosity_overshoot;
+
+        /*Collapsed formulation: */
+        Tria*  tria     =NULL;
+        Penta* pentabase=NULL;
+
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+        bool shelf;
+        Input* vx_input=NULL;
+        Input* vy_input=NULL;
+        Input* vxold_input=NULL;
+        Input* vyold_input=NULL;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
+
+        /*If on water, skip stiffness: */
+        if(onwater)return;
+
+        /*Find penta on bed as pattyn must be coupled to the dofs on the bed: */
+        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, numgrids);
+        GetDofList(&doflist);  //Pattyn dof list
+        GetDofList(&doflist2); //MacAyeal dof list
+
+        /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
+          get tria gaussian points as well as segment gaussian points. For tria gaussian 
+          points, order of integration is 2, because we need to integrate the product tB*D*B' 
+          which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
+          points, same deal, which yields 3 gaussian points.*/
+
+        order_area_gauss=5;
+        num_vert_gauss=5;
+
+        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
+
+        /*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: */
+        for (ig1=0; ig1<num_area_gauss; ig1++){
+                for (ig2=0; ig2<num_vert_gauss; ig2++){
+
+                        /*Pick up the gaussian point: */
+                        gauss_weight1=*(area_gauss_weights+ig1);
+                        gauss_weight2=*(vert_gauss_weights+ig2);
+                        gauss_weight=gauss_weight1*gauss_weight2;
+
+                        gauss_coord[0]=*(first_gauss_area_coord+ig1); 
+                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
+                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
+                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
+
+                        /*Get strain rate from velocity: */
+                        this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input);
+                        this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input);
+
+                        /*Get viscosity: */
+                        matice->GetViscosity3d(&viscosity, &epsilon[0]);
+                        matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
+
+                        /*Get B and Bprime matrices: */
+                        GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss_coord);
+                        tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
+
+                        /* Get Jacobian determinant: */
+                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
+
+                        /*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=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,numdof,1,
+                                                &D[0][0],3,3,0,
+                                                &Bprime[0][0],3,numdof2,0,
+                                                &Ke_gg_gaussian[0][0],0);
+
+                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
+                        for( i=0; i<numdof; i++){
+                                for(j=0;j<numdof2; j++){
+                                        Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+                                }
+                        }
+                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
+        } //for (ig1=0; ig1<num_area_gauss; ig1++)
+
+        /*Add Ke_gg to global matrix Kgg: */
+        // one need to be transposed
+        MatSetValues(Kgg,numdof,doflist,numdof2,doflist2,(const double*)Ke_gg,ADD_VALUES);
+        MatSetValues(Kgg,numdof2,doflist2,numdof,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->CreateKMatrixDiagnosticHorizFriction(Kgg);
+        }
+
+        delete tria;
+        delete pentabase;
+
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&fourth_gauss_vert_coord);
+        xfree((void**)&area_gauss_weights);
+        xfree((void**)&vert_gauss_weights);
+        xfree((void**)&first_gauss_area_coord2d);
+        xfree((void**)&second_gauss_area_coord2d);
+        xfree((void**)&third_gauss_area_coord2d);
+        xfree((void**)&gauss_weights2d);
+        xfree((void**)&doflist);
+}
+/*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/
 void  Penta::CreateKMatrixDiagnosticHutter(Mat Kgg){
@@ -2268 +2481 @@
         double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
-        double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
         double Jdet;
 
@@ -2798 +3010 @@
                                                 &Ke_gg_gaussian[0][0],0);
 
-                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
+                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
                         for( i=0; i<numdof; i++){
                                 for(j=0;j<numdof;j++){
@@ -3109 +3321 @@
                                         D_scalar=D_scalar*dt;
                                 }
-                                K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
-                                K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
+                                //K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
+                                //K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
+                                K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*(u)*(v);
+                                K[1][0]=D_scalar*(u)*(v);        K[1][1]=D_scalar*pow(v,2);
 
                                 /*Get B_artdiff: */

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

@@ -115 +115 @@
                 void      CreateKMatrixBalancedthickness(Mat Kggg);
                 void      CreateKMatrixBalancedvelocities(Mat Kggg);
+                void      CreateKMatrixCouplingMacAyealPattyn( Mat Kgg);
                 void      CreateKMatrixDiagnosticHutter( Mat Kgg);
                 void      CreateKMatrixDiagnosticMacAyeal( Mat Kgg);

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

@@ -52 +52 @@
 
 /*Reference Element numerics*/
+/*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/
+void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, double* gauss){
+        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
+         * For grid i, Bi can be expressed in the actual coordinate system
+         * by: 
+         *       Bi=[ dh/dx          0      ]
+         *          [   0           dh/dy   ]
+         *          [ 1/2*dh/dy  1/2*dh/dx  ]
+         * where h is the interpolation function for grid i.
+         *
+         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
+         */
+
+        int i;
+        const int numgrids=6;
+        const int NDOF3=3;
+        const int NDOF2=2;
+
+        double dh1dh6[NDOF3][numgrids];
+
+        /*Get dh1dh6 in actual coordinate system: */
+        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
+
+        /*Build B: */
+        for (i=0;i<numgrids;i++){
+                *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i]; 
+                *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
+
+                *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
+                *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
+
+                *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
+                *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 
+
+        }
+
+}
+/*}}}*/
 /*FUNCTION PentaRef::GetBPattyn {{{1*/
 void PentaRef::GetBPattyn(double* B, double* xyz_list, double* gauss){

issm/trunk/src/c/objects/Elements/PentaRef.h

@@ -32 +32 @@
                 void GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss);
                 void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
+                void GetBMacAyealPattyn(double* B, double* xyz_list, double* gauss);
                 void GetBPattyn(double* B, double* xyz_list, double* gauss);
                 void GetBprimePattyn(double* B, double* xyz_list, double* gauss);