Changeset 1047

Timestamp:

06/22/09 14:10:23 (16 years ago)

Author:

Mathieu Morlighem

Message:

cosmetics (bad commit)

File:

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

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

@@ -4 +4 @@
 
 #ifdef HAVE_CONFIG_H
-#include "config.h"
+        #include "config.h"
 #else
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
@@ -33 +33 @@
 
 Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],double tria_melting[3],
-                        double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,
-                        double tria_viscosity_overshoot,int tria_artdiff){
-
+                                double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,
+                                double tria_viscosity_overshoot,int tria_artdiff){
+        
         int i;
-
+        
         id=tria_id;
         mid=tria_mid;
@@ -66 +66 @@
         viscosity_overshoot=tria_viscosity_overshoot;
         artdiff=tria_artdiff;
-
+        
         return;
 }
@@ -156 +156 @@
         /*get enum type of Tria: */
         enum_type=TriaEnum();
-
+        
         /*marshall enum: */
         memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
-
+        
         /*marshall Tria data: */
         memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
@@ -187 +187 @@
         memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
         memcpy(marshalled_dataset,&artdiff,sizeof(artdiff));marshalled_dataset+=sizeof(artdiff);
-
+        
         *pmarshalled_dataset=marshalled_dataset;
         return;
 }
-
+                
 int   Tria::MarshallSize(){
         return sizeof(id)
-          +sizeof(mid)
-          +sizeof(mparid)
-          +sizeof(node_ids)
-          +sizeof(nodes)
-          +sizeof(node_offsets)
-          +sizeof(matice)
-          +sizeof(matice_offset)
-          +sizeof(matpar)
-          +sizeof(matpar_offset)
-          +sizeof(h)
-          +sizeof(s)
-          +sizeof(b)
-          +sizeof(k)
-          +sizeof(melting)
-          +sizeof(accumulation)
-          +sizeof(geothermalflux)
-          +sizeof(friction_type)
-          +sizeof(onbed)
-          +sizeof(p)
-          +sizeof(q)
-          +sizeof(shelf)
-          +sizeof(meanvel)
-          +sizeof(epsvel)
-          +sizeof(viscosity_overshoot)
-          +sizeof(artdiff)
-          +sizeof(int); //sizeof(int) for enum type
+                +sizeof(mid)
+                +sizeof(mparid)
+                +sizeof(node_ids)
+                +sizeof(nodes)
+                +sizeof(node_offsets)
+                +sizeof(matice)
+                +sizeof(matice_offset)
+                +sizeof(matpar)
+                +sizeof(matpar_offset)
+                +sizeof(h)
+                +sizeof(s)
+                +sizeof(b)
+                +sizeof(k)
+                +sizeof(melting)
+                +sizeof(accumulation)
+                +sizeof(geothermalflux)
+                +sizeof(friction_type)
+                +sizeof(onbed)
+                +sizeof(p)
+                +sizeof(q)
+                +sizeof(shelf)
+                +sizeof(meanvel)
+                +sizeof(epsvel)
+                +sizeof(viscosity_overshoot)
+                +sizeof(artdiff)
+                +sizeof(int); //sizeof(int) for enum type
 }
 
@@ -291 +291 @@
 
         int i;
-
+        
         DataSet* loadsin=NULL;
         DataSet* nodesin=NULL;
@@ -303 +303 @@
         /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */
         ResolvePointers((Object**)nodes,node_ids,node_offsets,3,nodesin);
-
+        
         /*Same for materials: */
         ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin);
@@ -317 +317 @@
         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
         if (analysis_type==ControlAnalysisEnum()){
-
+                
                 CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==DiagnosticAnalysisEnum()){
-
+        
                 if (sub_analysis_type==HorizAnalysisEnum()){
 
@@ -361 +361 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-
+        
         /* gaussian points: */
         int     num_gauss,ig;
@@ -375 +375 @@
         double newviscosity; //viscosity
         double oldviscosity; //viscosity
-
+        
         /* strain rate: */
         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
@@ -389 +389 @@
         double Ke_gg[numdof][numdof]; //local element stiffness matrix 
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
-
+        
         double Jdet;
-
+        
         /*input parameters for structural analysis (diagnostic): */
         double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
@@ -414 +414 @@
         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("El id %i Rank %i TriaElemnet input list before gaussian loop: \n",ELID,RANK); 
                 printf("   rho_ice: %g \n",matpar->GetRhoIce());
                 printf("   gravity: %g \n",matpar->GetG())
-                  printf("   rho_water: %g \n",matpar->GetRhoWater());
+                printf("   rho_water: %g \n",matpar->GetRhoWater());
                 printf("   Velocity: \n");
                 for (i=0;i<numgrids;i++){
@@ -430 +430 @@
                 printf("   bed [%g %g %g]\n",b[0],b[1],b[2]);
         }
-#endif
+        #endif
 
 
@@ -436 +436 @@
         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("   gaussian points: \n");
@@ -443 +443 @@
                 }
         }
-#endif
+        #endif
 
         /* Start  looping on the number of gaussian points: */
@@ -460 +460 @@
                 GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
                 GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
-
+                
                 /*Get viscosity: */
                 matice->GetViscosity2d(&viscosity, &epsilon[0]);
                 matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
-
+                
                 /* Get Jacobian determinant: */
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
 
                 /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
-                        onto this scalar matrix, so that we win some computational time: */
+                   onto this scalar matrix, so that we win some computational time: */
                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
                 D_scalar=newviscosity*thickness*gauss_weight*Jdet;
@@ -476 +476 @@
                         D[i][i]=D_scalar;
                 }
-
-#ifdef _DEBUGELEMENTS_
+                
+                #ifdef _DEBUGELEMENTS_
                 if(my_rank==RANK && id==ELID){ 
                         printf("   gaussian loop %i\n",ig);
@@ -492 +492 @@
                         printf("      gauss_weight: %g \n",gauss_weight);
                 }
-#endif
+                #endif
 
                 /*Get B and Bprime matrices: */
@@ -500 +500 @@
                 /*  Do the triple product tB*D*Bprime: */
                 TripleMultiply( &B[0][0],3,numdof,1,
-                                        &D[0][0],3,3,0,
-                                        &Bprime[0][0],3,numdof,0,
-                                        &Ke_gg_gaussian[0][0],0);
+                                          &D[0][0],3,3,0,
+                                          &Bprime[0][0],3,numdof,0,
+                                          &Ke_gg_gaussian[0][0],0);
 
                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
-
-#ifdef _DEBUGELEMENTS_
+                
+                #ifdef _DEBUGELEMENTS_
                 if(my_rank==RANK && id==ELID){ 
                         printf("      B:\n");
@@ -524 +524 @@
                         }
                 }
-#endif
+                #endif
         } // for (ig=0; ig<num_gauss; ig++)
 
@@ -537 +537 @@
 
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("      Ke_gg erms:\n");
@@ -552 +552 @@
 
         }
-#endif
-
-cleanup_and_return: 
+        #endif
+
+        cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -576 +576 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-
+        
         /* gaussian points: */
         int     num_gauss,ig;
@@ -597 +597 @@
         double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
         double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
-
+        
         double Jdettria;
-
+        
         /*input parameters for structural analysis (diagnostic): */
         double  vxvy_list[numgrids][2]={0.0};
@@ -628 +628 @@
                 vy_list[i]=vxvy_list[i][1];
         }
-
+        
         found=inputs->Recover("dt",&dt);
         if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
@@ -645 +645 @@
                 v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
                 v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
-
+                
                 K[0][0]=pow(Jdettria,.5)/2.0*fabs(v_gauss[0]);
                 K[1][1]=pow(Jdettria,.5)/2.0*fabs(v_gauss[1]);
@@ -671 +671 @@
                 /*  Do the triple product tL*D*L: */
                 TripleMultiply( &L[0],1,numdof,1,
-                                        &DL_scalar,1,1,0,
-                                        &L[0],1,numdof,0,
-                                        &Ke_gg_gaussian[0][0],0);
-
+                                          &DL_scalar,1,1,0,
+                                          &L[0],1,numdof,0,
+                                          &Ke_gg_gaussian[0][0],0);
+                
                 /*Get B  and B prime matrix: */
                 GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
                 GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
-
+                
                 //Get vx, vy and their derivatives at gauss point
                 GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
                 GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
-
+                
                 GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
                 GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
@@ -702 +702 @@
 
                 TripleMultiply( &B[0][0],2,numdof,1,
-                                        &DL[0][0],2,2,0,
-                                        &B[0][0],2,numdof,0,
-                                        &Ke_gg_thickness1[0][0],0);
+                                          &DL[0][0],2,2,0,
+                                          &B[0][0],2,numdof,0,
+                                          &Ke_gg_thickness1[0][0],0);
 
                 TripleMultiply( &B[0][0],2,numdof,1,
-                                        &DLprime[0][0],2,2,0,
-                                        &Bprime[0][0],2,numdof,0,
-                                        &Ke_gg_thickness2[0][0],0);
+                                          &DLprime[0][0],2,2,0,
+                                          &Bprime[0][0],2,numdof,0,
+                                          &Ke_gg_thickness2[0][0],0);
 
                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
@@ -715 +715 @@
                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
-
+                
                 if(artdiff){
-
+                        
                         /* Compute artificial diffusivity */
                         KDL[0][0]=DL_scalar*K[0][0];
@@ -723 +723 @@
 
                         TripleMultiply( &Bprime[0][0],2,numdof,1,
-                                                &KDL[0][0],2,2,0,
-                                                &Bprime[0][0],2,numdof,0,
-                                                &Ke_gg_gaussian[0][0],0);
+                                                  &KDL[0][0],2,2,0,
+                                                  &Bprime[0][0],2,numdof,0,
+                                                  &Ke_gg_gaussian[0][0],0);
 
                         /* Add artificial diffusivity matrix */
                         for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
-
-                }
-
-#ifdef _DEBUGELEMENTS_
+                        
+                }
+
+                #ifdef _DEBUGELEMENTS_
                 if(my_rank==RANK && id==ELID){ 
                         printf("      B:\n");
@@ -749 +749 @@
                         }
                 }
-#endif
+                #endif
         } // for (ig=0; ig<num_gauss; ig++)
 
@@ -755 +755 @@
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("      Ke_gg erms:\n");
@@ -770 +770 @@
 
         }
-#endif
-
-cleanup_and_return: 
+        #endif
+
+        cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -795 +795 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-
+        
         /* gaussian points: */
         int     num_gauss,ig;
@@ -809 +809 @@
         double L[numgrids];
         double Jdettria;
-
+        
         /*input parameters for structural analysis (diagnostic): */
         double  accumulation_list[numgrids]={0.0};
@@ -861 +861 @@
                 GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
                 GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3);
-
+                
                 /* Add value into pe_g: */
                 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
-
+                
         } // for (ig=0; ig<num_gauss; ig++)
 
@@ -870 +870 @@
         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
-cleanup_and_return: 
+        cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -893 +893 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-
+        
         /* gaussian points: */
         int     num_gauss,ig;
@@ -911 +911 @@
         double Ke_gg[numdof][numdof]; //local element stiffness matrix 
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
-
+        
         double Jdet;
-
+        
         /*slope: */
         double  slope[2]={0.0,0.0};
@@ -933 +933 @@
         /*recover pointers: */
         inputs=(ParameterInputs*)vinputs;
-
+        
         /* Get node coordinates and dof list: */
         GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
@@ -953 +953 @@
         /*Build alpha2_list used by drag stiffness matrix*/
         Friction* friction=NewFriction();
-
+        
         /*Initialize all fields: */
         friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
         strcpy(friction->element_type,"2d");
-
+        
         friction->gravity=matpar->GetG();
         friction->rho_ice=matpar->GetRhoIce();
@@ -974 +974 @@
         DeleteFriction(&friction);
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
         }
-#endif
+        #endif
 
         /* 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);
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("   gaussian points: \n");
@@ -990 +990 @@
                 }
         }
-#endif
+        #endif
 
         /* Start  looping on the number of gaussian points: */
@@ -1025 +1025 @@
                         DL[i][i]=DL_scalar;
                 }
-
+                
                 /*  Do the triple producte tL*D*L: */
                 TripleMultiply( &L[0][0],2,numdof,1,
@@ -1039 +1039 @@
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
-cleanup_and_return: 
+        cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -1062 +1062 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-
+        
         /* gaussian points: */
         int     num_gauss,ig;
@@ -1079 +1079 @@
         double Ke_gg[numdof][numdof]; //local element stiffness matrix 
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
-
+        
         double Jdet;
-
+        
         /* Get node coordinates and dof list: */
         GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
@@ -1100 +1100 @@
                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
 
-
+                
                 /*Get L matrix: */
                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
@@ -1106 +1106 @@
                 /* Get Jacobian determinant: */
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
+                
                 DL_scalar=gauss_weight*Jdet;
 
@@ -1121 +1121 @@
         /*Add Ke_gg to global matrix Kgg: */
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
-
-cleanup_and_return: 
+                
+        cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -1132 +1132 @@
 #define __FUNCT__ "Tria::CreatePVector"
 void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
-
+        
         /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
         if (analysis_type==ControlAnalysisEnum()){
-
+                
                 CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
-
+        
         }
         else if (analysis_type==DiagnosticAnalysisEnum()){
                 if (sub_analysis_type==HorizAnalysisEnum()){
-
+                
                         CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
-
+                
                 }
                 else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
         }
         else if (analysis_type==SlopeComputeAnalysisEnum()){
-
+                
                 CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
         }
@@ -1174 +1174 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-
+        
         /* parameters: */
         double  plastic_stress; 
@@ -1215 +1215 @@
 
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("gravity %g\n",matpar->GetG());
@@ -1224 +1224 @@
                 printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]);
         }
-#endif
+        #endif
 
 
@@ -1230 +1230 @@
         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("   gaussian points: \n");
@@ -1237 +1237 @@
                 }
         }
-#endif
+        #endif
 
 
@@ -1251 +1251 @@
                 /*Compute thickness at gaussian point: */
                 GetParameterValue(&thickness, &h[0],gauss_l1l2l3);
-
+        
                 GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
-
+                
                 /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the 
                  * element itself: */
@@ -1262 +1262 @@
                 /* Get Jacobian determinant: */
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
-                /*Get nodal functions: */
+                
+                 /*Get nodal functions: */
                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
 
@@ -1270 +1270 @@
 
 
-#ifdef _DEBUGELEMENTS_
+                #ifdef _DEBUGELEMENTS_
                 if(my_rank==RANK && id==ELID){ 
                         printf("      gaussian %i\n",ig);
@@ -1280 +1280 @@
                         if(friction_type==1)printf("      plastic_stress(%g)\n",plastic_stress);
                 }
-#endif
+                #endif
 
                 /*Build pe_g_gaussian vector: */
@@ -1303 +1303 @@
         } //for (ig=0; ig<num_gauss; ig++)
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("      pe_g->terms\n",ig);
@@ -1315 +1315 @@
                 }
         }
-#endif
+        #endif
 
         /*Add pe_g to global vector pg: */
         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
-cleanup_and_return: 
+        cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -1342 +1342 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-
+        
         /* gaussian points: */
         int     num_gauss,ig;
@@ -1391 +1391 @@
 
                 GetParameterDerivativeValue(&slope[0], &param[0],&xyz_list[0][0], gauss_l1l2l3);
-
+                
                 /* Get Jacobian determinant: */
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
-                /*Get nodal functions: */
+                
+                 /*Get nodal functions: */
                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
 
@@ -1414 +1414 @@
         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
-cleanup_and_return: 
+        cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -1445 +1445 @@
         inputs->Recover("accumulation",&accumulation[0],1,dofs,3,(void**)nodes);
         inputs->Recover("geothermalflux",&geothermalflux[0],1,dofs,3,(void**)nodes);
-
+        
         //Update material if necessary
         if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){
@@ -1452 +1452 @@
                 matice->SetB(B_average);
         }
-
+        
         if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){
                 B_average=(B_list[0]+B_list[1]+B_list[2])/3.0;
@@ -1459 +1459 @@
 
 }
-
+                
 void  Tria::GetDofList(int* doflist,int* pnumberofdofspernode){
 
@@ -1465 +1465 @@
         int doflist_per_node[MAXDOFSPERNODE];
         int numberofdofspernode;
-
+        
         for(i=0;i<3;i++){
                 nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
@@ -1491 +1491 @@
 #define __FUNCT__ "Tria::GetParameterValue"
 void Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){
-
+        
         /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian 
          * point specifie by gauss_l1l2l3: */
-
+        
         /*nodal functions: */
         double l1l2l3[3];
@@ -1513 +1513 @@
 #define __FUNCT__ "Tria::GetParameterDerivativeValue"
 void Tria::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3){
-
+         
         const int NDOF2=2;
         const int numgrids=3;
@@ -1523 +1523 @@
          * p is a vector of size 2x1 already allocated.
          */
-
+        
         double dh1dh2dh3_basic[NDOF2][numgrids]; //nodal derivative functions in basic coordinate system.
 
@@ -1548 +1548 @@
         GetB(&B[0][0], xyz_list, gauss_l1l2l3);
 
-#ifdef _DEBUG_
+        #ifdef _DEBUG_
         printf("B for grid1 : [ %lf   %lf  \n",B[0][0],B[0][1]);
         printf("              [ %lf   %lf  \n",B[1][0],B[1][1]);
@@ -1560 +1560 @@
         printf("              [ %lf   %lf  \n",B[1][4],B[1][5]);
         printf("              [ %lf   %lf  ]\n",B[2][4],B[2][5]);
-
+                
         for (i=0;i<numgrids;i++){
                 printf("Velocity for grid %i %lf %lf\n",i,*(vxvy_list+2*i+0),*(vxvy_list+2*i+1));
         }
-#endif
+        #endif
 
         /*Multiply B by velocity, to get strain rate: */
         MatrixMultiply( &B[0][0],3,NDOF2*numgrids,0,
-                                velocity,NDOF2*numgrids,1,0,
-                                epsilon,0);
+                                      velocity,NDOF2*numgrids,1,0,
+                                                  epsilon,0);
 
 }
@@ -1581 +1581 @@
 
         double x1,x2,x3,y1,y2,y3;
-
+        
         x1=*(xyz_list+3*0+0);
         y1=*(xyz_list+3*0+1);
@@ -1596 +1596 @@
                 printf("%s%s\n",__FUNCT__," error message: negative jacobian determinant!");
         }
-
+        
 }
 
@@ -1607 +1607 @@
 
         double x1,x2,x3,y1,y2,y3,z1,z2,z3;
-
+        
         x1=*(xyz_list+3*0+0);
         y1=*(xyz_list+3*0+1);
@@ -1625 +1625 @@
                 printf("%s%s\n",__FUNCT__," error message: negative jacobian determinant!");
         }
-
+        
 }
 
@@ -1643 +1643 @@
          * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)
          */
-
+        
         int i;
         const int NDOF2=2;
@@ -1654 +1654 @@
         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3);
 
-#ifdef _DEBUG_ 
+        #ifdef _DEBUG_ 
         for (i=0;i<3;i++){
                 printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh2dh3_basic[0][i],dh1dh2dh3_basic[1][i]);
         }
-#endif
+        #endif
 
         /*Build B: */
@@ -1686 +1686 @@
          * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
          */
-
+        
         int i;
         const int NDOF2=2;
@@ -1726 +1726 @@
          * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
          */
-
+        
         int i;
         const int NDOF2=2;
@@ -1737 +1737 @@
         GetNodalFunctions(l1l2l3, gauss_l1l2l3);
 
-#ifdef _DELUG_ 
+        #ifdef _DELUG_ 
         for (i=0;i<3;i++){
                 printf("Node %i  h=%lf \n",i,l1l2l3[i]);
         }
-#endif
+        #endif
 
         /*Build L: */
@@ -1773 +1773 @@
          * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
          */
-
+        
         int i;
         const int NDOF1=1;
@@ -1784 +1784 @@
         GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
 
-#ifdef _DEBUG_ 
+        #ifdef _DEBUG_ 
         for (i=0;i<3;i++){
                 printf("Node %i  h=%lf \n",i,l1l2l3[i]);
         }
-#endif
+        #endif
 
         /*Build B_prog: */
@@ -1811 +1811 @@
          * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
          */
-
+        
         int i;
         const int NDOF1=1;
@@ -1834 +1834 @@
 #define __FUNCT__ "Tria::GetNodalFunctions"
 void Tria::GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3){
-
+        
         /*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -1851 +1851 @@
 #define __FUNCT__ "Tria::GetNodalFunctionsDerivativesBasic"
 void Tria::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3){
-
+        
         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
          * basic coordinate system: */
@@ -1885 +1885 @@
 #define __FUNCT__ "Tria::GetNodalFunctionsDerivativesParams"
 void Tria::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3,double* gauss_l1l2l3){
-
+        
         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
          * natural coordinate system) at the gaussian point. */
@@ -1915 +1915 @@
         const int NDOF2=2;
         const int numgrids=3;
-
+        
         /*Call Jacobian routine to get the jacobian:*/
         GetJacobian(Jinv, xyz_list, gauss_l1l2l3);
@@ -1921 +1921 @@
         /*Invert Jacobian matrix: */
         MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
-
+                
 }
 
@@ -1935 +1935 @@
         double x1,y1,x2,y2,x3,y3;
         double sqrt3=sqrt(3.0);
-
+        
         x1=*(xyz_list+numgrids*0+0);
         y1=*(xyz_list+numgrids*0+1);
@@ -1958 +1958 @@
 }
 
-
+                
 void  Tria::GetNodes(void** vpnodes){
         int i;
@@ -1967 +1967 @@
         }
 }
-
+                
 
 int Tria::GetOnBed(){
@@ -1979 +1979 @@
 }
 void          Tria::GetBedList(double* bed_list){
-
+        
         int i;
         for(i=0;i<3;i++)bed_list[i]=b[i];
 
 }
-
+  
 
 Object* Tria::copy() {
-
+        
         return new Tria(*this); 
 
@@ -1998 +1998 @@
 
         int i;
-
+        
         /* node data: */
         const int    numgrids=3;
@@ -2104 +2104 @@
                 throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit));
         }
-
+        
         /* 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);
 
-#ifdef _DEBUGELEMENTS_
+        #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
                 printf("   gaussian points: \n");
@@ -2115 +2115 @@
                 }
         }
-#endif
-
+        #endif
+        
         /* Start  looping on the number of gaussian points: */
         for (ig=0; ig<num_gauss; ig++){
@@ -2127 +2127 @@
                 /* Get Jacobian determinant: */
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-#ifdef _DEBUG_ 
+                #ifdef _DEBUG_ 
                 printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
-#endif
-
+                #endif
+                
                 /* Get nodal functions value at gaussian point:*/
                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
@@ -2178 +2178 @@
                         throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit));
                 }
-
+                
                 /*Add due_g_gaussian vector to due_g: */
                 for( i=0; i<numdof; i++){
@@ -2184 +2184 @@
                 }
         }
-
+        
         /*Add due_g to global vector du_g: */
         VecSetValues(du_g,numdof,doflist,(const double*)due_g,ADD_VALUES);
-
-cleanup_and_return: 
+        
+        cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -2252 +2252 @@
         double  lambda,mu;
         double  bed,thickness,Neff;
-
+        
         /*drag: */
         double  pcoeff,qcoeff;