Changeset 5847

Timestamp:

09/16/10 13:11:42 (15 years ago)

Author:

Mathieu Morlighem

Message:

split MacAyeal2d and MacAyeal3d in Penta

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (4 diffs)
• Penta.h (modified) (1 diff)

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

@@ -2193 +2193 @@
         switch(approximation){
                 case MacAyealApproximationEnum:
-                        CreateKMatrixDiagnosticMacAyeal( Kgg);
+                        CreateKMatrixDiagnosticMacAyeal2d( Kgg);
                         break;
                 case PattynApproximationEnum:
@@ -2206 +2206 @@
                         break;
                 case MacAyealPattynApproximationEnum:
-                        CreateKMatrixDiagnosticMacAyeal( Kgg);
+                        CreateKMatrixDiagnosticMacAyeal3d( Kgg);
                         CreateKMatrixDiagnosticPattyn( Kgg);
                         CreateKMatrixCouplingMacAyealPattyn( Kgg);
@@ -2290 +2290 @@
 
 }
-/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal{{{1*/
-void Penta::CreateKMatrixDiagnosticMacAyeal( Mat Kgg){
+/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal2d{{{1*/
+void Penta::CreateKMatrixDiagnosticMacAyeal2d(Mat Kgg){
         
-        /* local declarations */
-        int             i,j;
-
-        /* node data: */
+        /*If on water, skip stiffness: */
+        if(IsOnWater()) 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 (!IsOnBed()){
+                /*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{
+
+                /*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=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->CreateKMatrix(Kgg,NULL,NULL);
+                delete tria->matice; delete tria;
+
+                /*Delete B averaged*/
+                this->matice->inputs->DeleteInput(RheologyBbarEnum);
+
+                return;
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal3d{{{1*/
+void Penta::CreateKMatrixDiagnosticMacAyeal3d( Mat Kgg){
+
         const int    numdof2d=2*NUMVERTICES2D;
+        int             i,j,ig;
         double       xyz_list[NUMVERTICES][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 viscosity, oldviscosity, newviscosity, viscosity_overshoot;
         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];
@@ -2323 +2344 @@
         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;
-        int  approximation;
-
-        inputs->GetParameterValue(&approximation,ApproximationEnum);
 
         /*If on water, skip stiffness: */
-        if(approximation==MacAyealApproximationEnum){
-                if(IsOnWater())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 (!IsOnBed()){
-                        /*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{
-
-                        /*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,NULL,NULL);
-                        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(IsOnWater())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,NUMVERTICES);
-                tria->GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum);  //Pattyn dof list
-
-                /*Retrieve all inputs we will be needing: */
-                Input* vx_input=inputs->GetInput(VxEnum);       ISSMASSERT(vx_input);
-                Input* vy_input=inputs->GetInput(VyEnum);       ISSMASSERT(vy_input);
-                Input* vxold_input=inputs->GetInput(VxOldEnum); ISSMASSERT(vxold_input);
-                Input* vyold_input=inputs->GetInput(VyOldEnum); ISSMASSERT(vyold_input);
-
-                /* 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(IsOnBed() && !IsOnShelf()){
-                        /*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,NULL,NULL);
-                }
-
-                /*Clean up and return*/
-                delete tria->matice;
-                delete tria;
-                delete gauss_tria;
-                delete gauss;
-                xfree((void**)&doflist);
-        }
+        if(IsOnWater())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 dof list: */
+        tria->GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum);  //Pattyn dof list
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES);
+        this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
+        Input* vx_input=inputs->GetInput(VxEnum);       ISSMASSERT(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);       ISSMASSERT(vy_input);
+        Input* vxold_input=inputs->GetInput(VxOldEnum); ISSMASSERT(vxold_input);
+        Input* vyold_input=inputs->GetInput(VyOldEnum); ISSMASSERT(vyold_input);
+
+        /* 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);
+
+                tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria);
+                tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
+
+
+                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);
+                matice->GetViscosity3d(&viscosity, &epsilon[0]);
+                matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
+                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                D_scalar=2*newviscosity*gauss->weight*Jdet;
+                for (i=0;i<3;i++) D[i][i]=D_scalar;
+
+                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);
+
+                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+        }
+
+        MatSetValues(Kgg,numdof2d,doflist,numdof2d,doflist,(const double*)Ke_gg,ADD_VALUES);
+
+        //Deal with 2d friction at the bedrock interface
+        if(IsOnBed() && !IsOnShelf()){
+                /*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,NULL,NULL);
+        }
+
+        /*Clean up and return*/
+        delete tria->matice;
+        delete tria;
+        delete gauss_tria;
+        delete gauss;
+        xfree((void**)&doflist);
 }
 /*}}}*/

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

@@ -127 +127 @@
                 void      CreateKMatrixDiagnosticHoriz( Mat Kgg);
                 void      CreateKMatrixDiagnosticHutter( Mat Kgg);
-                void      CreateKMatrixDiagnosticMacAyeal( Mat Kgg);
+                void      CreateKMatrixDiagnosticMacAyeal2d(Mat Kgg);
+                void      CreateKMatrixDiagnosticMacAyeal3d(Mat Kgg);
                 void      CreateKMatrixDiagnosticPattyn( Mat Kgg);
                 void      CreateKMatrixDiagnosticStokes( Mat Kgg);