Changeset 5387 for issm/trunk

Timestamp:

08/18/10 16:08:34 (15 years ago)

Author:

seroussi

Message:

moved HorizFriction into MacAyeal, Pattyn and MacAyealPattynCOupling friction

Location:

issm/trunk/src/c/objects/Elements

Files:

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

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

@@ -2502 +2502 @@
                  * grids: */
 
-                tria->CreateKMatrixDiagnosticHorizFriction(Kgg);
+                tria->CreateKMatrixCouplingMacAyealPattynFriction(Kgg);
         }
 
@@ -2835 +2835 @@
 
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreateKMatrixDiagnosticHorizFriction(Kgg);
+                tria->CreateKMatrixDiagnosticPattynFriction(Kgg);
                 delete tria->matice; delete tria;
         }

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

@@ -3591 +3591 @@
         /*Do not forget to include friction: */
         if(!shelf){
-                CreateKMatrixDiagnosticHorizFriction(Kgg);
+                CreateKMatrixDiagnosticMacAyealFriction(Kgg);
         }
 
@@ -3601 +3601 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixDiagnosticHorizFriction {{{1*/
-void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg){
+/*FUNCTION Tria::CreateKMatrixCouplingMacAyealPattynFriction {{{1*/
+void  Tria::CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg){
 
         /* local declarations */
@@ -3612 +3612 @@
         const int numdof   = 2 *numvertices;
         double    xyz_list[numvertices][3];
-        int*      doflist=NULL;
+        int*      doflistm=NULL;
+        int*      doflistp=NULL;
         int       numberofdofspernode=2;
 
@@ -3669 +3670 @@
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
-        GetDofList(&doflist);
+        GetDofList(&doflistm,MacAyealApproximationEnum);
+        GetDofList(&doflistp,PattynApproximationEnum);
+
+        if (shelf){
+                /*no friction, do nothing*/
+                return;
+        }
+
+        /*build friction object, used later on: */
+        if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
+        friction=new Friction("2d",inputs,matpar,analysis_type);
+
+        /*COMPUT alpha2_list (TO BE DELETED)*/
+        for(i=0;i<numvertices;i++){
+                friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum);
+        }
+
+        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+        /* Start  looping on the number of gaussian points: */
+        for (ig=0; ig<num_gauss; ig++){
+                /*Pick up the gaussian point: */
+                gauss_weight=*(gauss_weights+ig);
+                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+                /*Friction: */
+                // friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
+                TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED
+
+                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
+                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
+                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
+                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+
+                if (slope_magnitude>MAXSLOPE){
+                        alpha2=pow((double)10,MOUNTAINKEXPONENT);
+                }
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+
+                /*Get L matrix: */
+                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+                
+                DL_scalar=alpha2*gauss_weight*Jdet;
+                for (i=0;i<2;i++){
+                        DL[i][i]=DL_scalar;
+                }
+                
+                /*  Do the triple producte tL*D*L: */
+                TripleMultiply( &L[0][0],2,numdof,1,
+                                        &DL[0][0],2,2,0,
+                                        &L[0][0],2,numdof,0,
+                                        &Ke_gg_gaussian[0][0],0);
+
+                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+
+        } // for (ig=0; ig<num_gauss; ig++)
+
+        /*Add Ke_gg to global matrix Kgg: */
+        MatSetValues(Kgg,numdof,doflistm,numdof,doflistp,(const double*)Ke_gg,ADD_VALUES);
+        MatSetValues(Kgg,numdof,doflistp,numdof,doflistm,(const double*)Ke_gg,ADD_VALUES);
+
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+        xfree((void**)&doflistm);
+        xfree((void**)&doflistp);
+        delete friction;
+}       
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
+void  Tria::CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg){
+
+        /* local declarations */
+        int       i,j;
+        int       analysis_type;
+
+        /* node data: */
+        const int numvertices = 3;
+        const int numdof   = 2 *numvertices;
+        double    xyz_list[numvertices][3];
+        int*      doflist=NULL;
+        int       numberofdofspernode=2;
+
+        /* 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 *gauss_weights           = NULL;
+        double  gauss_weight;
+        double  gauss_l1l2l3[3];
+
+        /* matrices: */
+        double L[2][numdof];
+        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
+        double DL_scalar;
+
+        /* local element matrices: */
+        double Ke_gg[numdof][numdof]={0.0};
+        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
+        
+        double Jdet;
+        
+        /*slope: */
+        double  slope[2]={0.0,0.0};
+        double  slope_magnitude;
+
+        /*friction: */
+        Friction *friction = NULL;
+        double    alpha2;
+        double    alpha2_list[numvertices];                                       //TO BE DELETED
+        double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
+
+        double MAXSLOPE=.06; // 6 %
+        double MOUNTAINKEXPONENT=10;
+
+        /*inputs: */
+        bool shelf;
+        int  drag_type;
+        Input* surface_input=NULL;
+        Input* vx_input=NULL;
+        Input* vy_input=NULL;
+        Input* vz_input=NULL;
+
+        /*retrive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+        inputs->GetParameterValue(&drag_type,DragTypeEnum);
+        surface_input=inputs->GetInput(SurfaceEnum);
+        vx_input=inputs->GetInput(VxEnum);
+        vy_input=inputs->GetInput(VyEnum);
+        vz_input=inputs->GetInput(VzEnum);
+        
+        /* Get node coordinates and dof list: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+        GetDofList(&doflist,MacAyealApproximationEnum);
+
+        if (shelf){
+                /*no friction, do nothing*/
+                return;
+        }
+
+        /*build friction object, used later on: */
+        if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
+        friction=new Friction("2d",inputs,matpar,analysis_type);
+
+        /*COMPUT alpha2_list (TO BE DELETED)*/
+        for(i=0;i<numvertices;i++){
+                friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum);
+        }
+
+        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+        /* Start  looping on the number of gaussian points: */
+        for (ig=0; ig<num_gauss; ig++){
+                /*Pick up the gaussian point: */
+                gauss_weight=*(gauss_weights+ig);
+                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+                /*Friction: */
+                // friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
+                TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED
+
+                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
+                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
+                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
+                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+
+                if (slope_magnitude>MAXSLOPE){
+                        alpha2=pow((double)10,MOUNTAINKEXPONENT);
+                }
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+
+                /*Get L matrix: */
+                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+                
+                DL_scalar=alpha2*gauss_weight*Jdet;
+                for (i=0;i<2;i++){
+                        DL[i][i]=DL_scalar;
+                }
+                
+                /*  Do the triple producte tL*D*L: */
+                TripleMultiply( &L[0][0],2,numdof,1,
+                                        &DL[0][0],2,2,0,
+                                        &L[0][0],2,numdof,0,
+                                        &Ke_gg_gaussian[0][0],0);
+
+                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+
+        } // for (ig=0; ig<num_gauss; ig++)
+
+        /*Add Ke_gg to global matrix Kgg: */
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
+
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+        xfree((void**)&doflist);
+        delete friction;
+}       
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixDiagnosticPattynFriction {{{1*/
+void  Tria::CreateKMatrixDiagnosticPattynFriction(Mat Kgg){
+
+        /* local declarations */
+        int       i,j;
+        int       analysis_type;
+
+        /* node data: */
+        const int numvertices = 3;
+        const int numdof   = 2 *numvertices;
+        double    xyz_list[numvertices][3];
+        int*      doflist=NULL;
+        int       numberofdofspernode=2;
+
+        /* 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 *gauss_weights           = NULL;
+        double  gauss_weight;
+        double  gauss_l1l2l3[3];
+
+        /* matrices: */
+        double L[2][numdof];
+        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
+        double DL_scalar;
+
+        /* local element matrices: */
+        double Ke_gg[numdof][numdof]={0.0};
+        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
+        
+        double Jdet;
+        
+        /*slope: */
+        double  slope[2]={0.0,0.0};
+        double  slope_magnitude;
+
+        /*friction: */
+        Friction *friction = NULL;
+        double    alpha2;
+        double    alpha2_list[numvertices];                                       //TO BE DELETED
+        double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
+
+        double MAXSLOPE=.06; // 6 %
+        double MOUNTAINKEXPONENT=10;
+
+        /*inputs: */
+        bool shelf;
+        int  drag_type;
+        Input* surface_input=NULL;
+        Input* vx_input=NULL;
+        Input* vy_input=NULL;
+        Input* vz_input=NULL;
+
+        /*retrive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+        inputs->GetParameterValue(&drag_type,DragTypeEnum);
+        surface_input=inputs->GetInput(SurfaceEnum);
+        vx_input=inputs->GetInput(VxEnum);
+        vy_input=inputs->GetInput(VyEnum);
+        vz_input=inputs->GetInput(VzEnum);
+        
+        /* Get node coordinates and dof list: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+        GetDofList(&doflist,PattynApproximationEnum);
 
         if (shelf){

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

@@ -124 +124 @@
                 void      CreateKMatrixBalancedvelocities(Mat Kgg);
                 void      CreateKMatrixDiagnosticMacAyeal(Mat Kgg);
-                void      CreateKMatrixDiagnosticHorizFriction(Mat Kgg);
+                void      CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg);
+                void      CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg);
+                void      CreateKMatrixDiagnosticPattynFriction(Mat Kgg);
                 void      CreateKMatrixDiagnosticHutter(Mat Kgg);
                 void      CreateKMatrixDiagnosticSurfaceVert(Mat Kgg);