Changeset 5910

Timestamp:

09/20/10 15:14:23 (15 years ago)

Author:

seroussi

Message:

clean tria

File:

• issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (12 diffs)

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

@@ -2405 +2405 @@
 ElementMatrix* Tria::CreateKMatrixBalancedthickness_CG(void){
 
+
         /*Constants*/
         const int    numdof=NDOF1*NUMVERTICES;
@@ -3399 +3400 @@
 ElementMatrix* Tria::CreateKMatrixThermal(void){
 
+        /*Constants*/
         const int    numdof=NDOF1*NUMVERTICES;
-        int i,j;
-        int     ig;
-        double       xyz_list[NUMVERTICES][3];
-        double mixed_layer_capacity;
-        double thermal_exchange_velocity;
-        double rho_water;
-        double rho_ice;
-        double heatcapacity;
+
+        /*Intermediaries */
+        int       i,j,ig;
+        double    mixed_layer_capacity,thermal_exchange_velocity;
+        double    rho_ice,rho_water,heatcapacity;
+        double    Jdet,dt;
+        double    xyz_list[NUMVERTICES][3];
+        double    l1l2l3[NUMVERTICES];
+        double    tl1l2l3D[3];
+        double    D_scalar;
+        double    Ke_gaussian[numdof][numdof]={0.0};
         GaussTria *gauss=NULL;
-        double  Jdet;
-        double  K_terms[numdof][numdof]={0.0};
-        double  Ke_gaussian[numdof][numdof]={0.0};
-        double  l1l2l3[NUMVERTICES];
-        double     tl1l2l3D[3];
-        double  D_scalar;
-        double dt;
 
         /*Initialize Element matrix and return if necessary*/
@@ -3436 +3434 @@
                 gauss->GaussPoint(ig);
                 
+                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
                 GetNodalFunctions(&l1l2l3[0], gauss);
                                 
-                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
                 D_scalar=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
                 if(dt) D_scalar=dt*D_scalar;
@@ -3664 +3662 @@
 void Tria::CreatePVectorDiagnosticMacAyeal( Vec pg, Vec pf){
 
-        int             i,j,ig;
-
-        /* node data: */
+        /*Constants*/
         const int    numdof=NDOF2*NUMVERTICES;
-        double       xyz_list[NUMVERTICES][3];
-        
-        /* parameters: */
-        double  plastic_stress; 
-        double  slope[2];
-        double  driving_stress_baseline;
-        GaussTria* gauss=NULL;
-
-        double Jdet;
-        double l1l2l3[3];
-        double  pe_g[numdof]={0.0};
-        double  pe_g_gaussian[numdof];
-        double  thickness;
-        int  drag_type;
-
-        /*intermediary: */
+
+        /*Intermediaries */
+        int            i,j,ig,drag_type;
+        double         plastic_stress,driving_stress_baseline,thickness;
+        double         Jdet;
+        double         xyz_list[NUMVERTICES][3];
+        double         slope[2];
+        double         l1l2l3[3];
+        double         pe_g[numdof]={0.0};
+        double         pe_g_gaussian[numdof];
+        GaussTria*     gauss=NULL;
         ElementVector* pe=NULL;
 
-        /*retrieve inputs :*/
+        /*First, if we are on water, return empty vector: */
+        if(IsOnWater()) return;
+
+        /*Retrieve all Inputs and parameters: */
         inputs->GetParameterValue(&drag_type,DragTypeEnum);
         Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 
@@ -3692 +3686 @@
         Input* drag_input=inputs->GetInput(DragCoefficientEnum);ISSMASSERT(drag_input);
 
-        /*First, if we are on water, return empty vector: */
-        if(IsOnWater()) return;
-
-        /* Get node coordinates and dof list: */
+        /*Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
@@ -3703 +3694 @@
 
                 gauss->GaussPoint(ig);
+
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+                GetNodalFunctions(l1l2l3, gauss);
 
                 thickness_input->GetParameterValue(&thickness,gauss);
@@ -3709 +3703 @@
                 /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the 
                  * element itself: */
-                if(drag_type==1){
-                        drag_input->GetParameterValue(&plastic_stress,gauss);
-                }
-
-                /* Get Jacobian determinant: */
-                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-                
-                 /*Get nodal functions: */
-                GetNodalFunctions(l1l2l3, gauss);
-
-                /*Compute driving stress: */
+                if(drag_type==1) drag_input->GetParameterValue(&plastic_stress,gauss);
+
                 driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
 
@@ -3738 +3723 @@
                 }
 
-                /*Add pe_g_gaussian vector to pe_g: */
                 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
-
-        }
-
-        /*Initialize element vector: */
+        }
+
         pe=NewElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-
-        /*Add pe_g values to pe element stifness load: */
         pe->AddValues(&pe_g[0]);
-
-        /*Add pe element load vector to global load vector: */
         pe->AddToGlobal(pg,pf);
 
         /*Free ressources:*/
         delete pe;
-
-        /*Free ressources:*/
         delete gauss;
 
@@ -3764 +3740 @@
 void Tria::CreatePVectorAdjointBalancedthickness(Vec p_g){
 
-        int i;
-
-        /* node data: */
+        /*Constants*/
         const int    numdof=1*NUMVERTICES;
-        double       xyz_list[NUMVERTICES][3];
-        int*         doflist=NULL;
-
-        /* gaussian points: */
-        int     ig;
+
+        /*Intermediaries */
+        int         i,ig;
+        int*        doflist=NULL;
+        double      Jdet;
+        double      thickness,thicknessobs,weight;
+        double      xyz_list[NUMVERTICES][3];
+        double      l1l2l3[3];
+        double      pe_g_gaussian[numdof];
+        double      pe_g[numdof]                 ={0.0};
         GaussTria* gauss=NULL;
-
-        /* parameters: */
-        double  thickness,thicknessobs,weight;
-
-        /*element vector : */
-        double  pe_g[numdof]={0.0};
-        double  pe_g_gaussian[numdof];
-
-        /* Jacobian: */
-        double Jdet;
-
-        /*nodal functions: */
-        double l1l2l3[3];
 
         /* Get node coordinates and dof list: */
@@ -3792 +3758 @@
         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 
-        /*Retrieve all inputs we will be needing: */
+        /*Retrieve all Inputs and parameters: */
         Input* thickness_input   =inputs->GetInput(ThicknessEnum);   ISSMASSERT(thickness_input);
         Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input);
@@ -3803 +3769 @@
                 gauss->GaussPoint(ig);
 
-                /* Get Jacobian determinant: */
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-
-                /* Get nodal functions value at gaussian point:*/
                 GetNodalFunctions(l1l2l3, gauss);
 
-                /*Get parameters at gauss point*/
                 thickness_input->GetParameterValue(&thickness, gauss);
                 thicknessobs_input->GetParameterValue(&thicknessobs, gauss);
@@ -3819 +3781 @@
 
                 /*Add pe_g_gaussian vector to pe_g: */
-                for( i=0; i<numdof; i++){
-                        pe_g[i]+=pe_g_gaussian[i];
-                }
-        }
-
-        /*Add pe_g to global vector p_g: */
+                for(i=0; i<numdof;i++) pe_g[i]+=pe_g_gaussian[i];
+        }
+
         VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);