Changeset 8418

Timestamp:

05/24/11 15:23:22 (14 years ago)

Author:

seroussi

Message:

last moved from tria to penta : CreateKMatrixCouplingMacAyealPattynFriction

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (1 diff)
• Tria.cpp (modified) (1 diff)
• Tria.h (modified) (1 diff)

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

@@ -777 +777 @@
 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynFriction(void){
 
+        /*Constants*/
+        const int numdof        = NDOF2 *NUMVERTICES;
+        const int numdoftotal   = NDOF4 *NUMVERTICES;
+        
+        /*Intermediaries */
+        int       i,j,ig,analysis_type,drag_type;
+        double    Jdet2d,slope_magnitude,alpha2;
+        double    xyz_list[NUMVERTICES][3];
+        double    xyz_list_tria[NUMVERTICES2D][3]={0.0};
+        double    slope[3]={0.0,0.0,0.0};
+        double    MAXSLOPE=.06; // 6 %
+        double    MOUNTAINKEXPONENT=10;
+        double    L[2][numdof];
+        double    DL[2][2]                  ={{ 0,0 },{0,0}}; //for basal drag
+        double    DL_scalar;
+        double    Ke_gg[numdof][numdof]     ={0.0};
+        double    Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
+        Friction  *friction = NULL;
+        GaussPenta *gauss=NULL;
+
         /*Initialize Element matrix and return if necessary*/
         if(IsOnShelf() || !IsOnBed()) return NULL;
-
-        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-        ElementMatrix* Ke=tria->CreateKMatrixCouplingMacAyealPattynFriction();
-        delete tria->matice; delete tria;
-
+        ElementMatrix* Ke1=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
+        ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+        delete Ke1; delete Ke2;
+
+        /*retrieve inputs :*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        inputs->GetParameterValue(&drag_type,DragTypeEnum);
+        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
+        Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
+        Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
+
+        /*build friction object, used later on: */
+        if (drag_type!=2)_error_(" non-viscous friction not supported yet!");
+        friction=new Friction("2d",inputs,matpar,analysis_type);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussPenta(0,1,2,2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /*Friction: */
+                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
+
+                // 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);
+                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+
+                if (slope_magnitude>MAXSLOPE){
+                        alpha2=pow((double)10,MOUNTAINKEXPONENT);
+                }
+
+                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
+                GetL(&L[0][0], gauss,NDOF2);
+
+                DL_scalar=alpha2*gauss->weight*Jdet2d;
+                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(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i][j];
+        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i][j];
+
+        /*Clean up and return*/
+        delete gauss;
+        delete friction;
         return Ke;
 }

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

@@ -904 +904 @@
         return Ke;
 }
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixCouplingMacAyealPattynFriction {{{1*/
-ElementMatrix* Tria::CreateKMatrixCouplingMacAyealPattynFriction(void){
-
-        /*Constants*/
-        const int numdof        = NDOF2 *NUMVERTICES;
-        const int numdoftotal   = NDOF4 *NUMVERTICES;
-        
-        /*Intermediaries */
-        int       i,j,ig,analysis_type,drag_type;
-        double    Jdet,slope_magnitude,alpha2;
-        double    xyz_list[NUMVERTICES][3];
-        double    slope[2]={0.0,0.0};
-        double    MAXSLOPE=.06; // 6 %
-        double    MOUNTAINKEXPONENT=10;
-        double    L[2][numdof];
-        double    DL[2][2]                  ={{ 0,0 },{0,0}}; //for basal drag
-        double    DL_scalar;
-        double    Ke_gg[numdof][numdof]     ={0.0};
-        double    Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
-        Friction  *friction = NULL;
-        GaussTria *gauss=NULL;
-
-        /*Initialize Element matrix and return if necessary*/
-        if(IsOnShelf()) return NULL;
-        ElementMatrix* Ke1=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-        ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
-        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
-        delete Ke1; delete Ke2;
-
-        /*retrieve inputs :*/
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        inputs->GetParameterValue(&drag_type,DragTypeEnum);
-        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
-        Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
-        Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
-        Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
-
-        /*build friction object, used later on: */
-        if (drag_type!=2)_error_(" non-viscous friction not supported yet!");
-        friction=new Friction("2d",inputs,matpar,analysis_type);
-
-        /* Start  looping on the number of gaussian points: */
-        gauss=new GaussTria(2);
-        for (ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                /*Friction: */
-                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
-
-                // 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);
-                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);
-
-                /*Get L matrix: */
-                GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
-
-                
-                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(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i][j];
-        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i][j];
-
-        /*Clean up and return*/
-        delete gauss;
-        delete friction;
-        return Ke;
-}       
 /*}}}*/
 /*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/

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

@@ -150 +150 @@
                 ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void);
                 ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void);
-                ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
                 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
                 ElementMatrix* CreateKMatrixMelting(void);