Changeset 8417

Timestamp:

05/24/11 15:09:01 (14 years ago)

Author:

seroussi

Message:

moved CreateKMatrix DiagnosticPattynFriction from Tria to Penta

Location:

issm/trunk/src/c/objects/Elements

Files:

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

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

@@ -1371 +1371 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){
 
+        /*Constants*/
+        const int numdof   = NDOF2*NUMVERTICES;
+        
+        /*Intermediaries */
+        int       i,j,ig;
+        int       analysis_type,drag_type;
+        double    xyz_list[NUMVERTICES][3];
+        double    xyz_list_tria[NUMVERTICES2D][3]={0.0};
+        double    slope_magnitude,alpha2,Jdet;
+        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_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;
 
-        /*Build a tria element using the 3 nodes of the base of the penta. Then use 
-         * the tria functionality to build a friction stiffness matrix on these 3
-         * nodes: */
-        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-        ElementMatrix* Ke=tria->CreateKMatrixDiagnosticPattynFriction();
-        delete tria->matice; delete tria;
-
-        /*clean-up and return*/
+        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+
+        /*Retrieve all inputs and parameters*/
+        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);
+
+                GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss);
+                GetL(&L[0][0], gauss,NDOF2);
+
+                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
+                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); 
+                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+
+                // 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: */
+                if (slope_magnitude>MAXSLOPE){
+                        alpha2=pow((double)10,MOUNTAINKEXPONENT);
+                }
+                
+                DL_scalar=alpha2*gauss->weight*Jdet;
+                for (i=0;i<2;i++) DL[i][i]=DL_scalar;
+                
+                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->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
+        }
+
+        /*Clean up and return*/
+        delete gauss;
+        delete friction;
         return Ke;
 }

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

@@ -519 +519 @@
 
 }
+/*}}}*/
+/*FUNCTION PentaRef::GetL{{{1*/
+void PentaRef::GetL(double* L, GaussPenta* gauss, int numdof){
+        /*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+         ** For node i, Li can be expressed in the actual coordinate system
+         ** by: 
+         **       numdof=1: 
+         **                 Li=h;
+         **       numdof=2:
+         **                 Li=[ h   0 ]
+         **                    [ 0   h ]
+         ** where h is the interpolation function for node i.
+         **
+         ** We assume L has been allocated already, of size: NUMNODESP1 (numdof=1), or numdofx(numdof*NUMNODESP1) (numdof=2)
+         **/
+
+        int i;
+        double l1l6[6];
+
+        /*Get l1l6 in actual coordinate system: */
+        GetNodalFunctionsP1(l1l6,gauss);
+
+        /*Build L: */
+        if(numdof==1){
+                for (i=0;i<NUMNODESP1;i++){
+                        L[i]=l1l6[i]; 
+                }
+        }
+        else{
+                for (i=0;i<NUMNODESP1;i++){
+                        *(L+numdof*NUMNODESP1*0+numdof*i)=l1l6[i]; 
+                        *(L+numdof*NUMNODESP1*0+numdof*i+1)=0;
+                        *(L+numdof*NUMNODESP1*1+numdof*i)=0;
+                        *(L+numdof*NUMNODESP1*1+numdof*i+1)=l1l6[i];
+                }
+        }
+} 
 /*}}}*/
 /*FUNCTION PentaRef::GetLStokes {{{1*/

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

@@ -49 +49 @@
                 void GetBVert(double* B, double* xyz_list, GaussPenta* gauss);
                 void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss);
+                void GetL(double* L, GaussPenta* gauss,int numdof);
                 void GetLStokes(double* LStokes, GaussPenta* gauss);
                 void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss);

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

@@ -988 +988 @@
         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::CreateKMatrixDiagnosticPattynFriction {{{1*/
-ElementMatrix* Tria::CreateKMatrixDiagnosticPattynFriction(void){
-
-        /*Constants*/
-        const int numdof   = NDOF2*NUMVERTICES;
-        
-        /*Intermediaries */
-        int       i,j,ig;
-        int       analysis_type,drag_type;
-        double    xyz_list[NUMVERTICES][3];
-        double    slope_magnitude,alpha2,Jdet;
-        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_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* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
-
-        /*Retrieve all inputs and parameters*/
-        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);
-
-                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-                GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
-
-                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
-                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
-                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
-
-                // 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: */
-                if (slope_magnitude>MAXSLOPE){
-                        alpha2=pow((double)10,MOUNTAINKEXPONENT);
-                }
-                
-                DL_scalar=alpha2*gauss->weight*Jdet;
-                for (i=0;i<2;i++) DL[i][i]=DL_scalar;
-                
-                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->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
-        }
 
         /*Clean up and return*/

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

@@ -151 +151 @@
                 ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void);
                 ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
-                ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
                 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
                 ElementMatrix* CreateKMatrixMelting(void);