Changeset 511

Timestamp:

05/20/09 07:46:55 (16 years ago)

Author:

Mathieu Morlighem

Message:

first fixing of thermal runs

Location:

issm/trunk/src/c/objects

Files:

• Penta.cpp (modified) (153 diffs)
• Tria.cpp (modified) (5 diffs)

issm/trunk/src/c/objects/Penta.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!"
@@ -20 +20 @@
 }
 Penta::Penta( int penta_id, int penta_mid, int penta_mparid, int penta_node_ids[6], double penta_h[6], double penta_s[6], double penta_b[6], double penta_k[6], int penta_friction_type, 
-                                double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel, 
-                                int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], 
-                                int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_stokesreconditioning){
-        
+                        double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel, 
+                        int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], 
+                        int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_stokesreconditioning){
+
         int i;
 
@@ -83 +83 @@
         printf("   b=[%i,%i,%i,%i,%i,%i]\n",b[0],b[1],b[2],b[3],b[4],b[5]);
         printf("   k=[%i,%i,%i,%i,%i,%i]\n",k[0],k[1],k[2],k[3],k[4],k[5]);
-        
+
         printf("   friction_type: %i\n",friction_type);
         printf("   p: %g\n",p);
@@ -93 +93 @@
         printf("   epsvel: %g\n",epsvel);
         printf("   collapse: %i\n",collapse);
-        
+
         printf("   melting=[%i,%i,%i,%i,%i,%i]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]);
         printf("   accumulation=[%i,%i,%i,%i,%i,%i]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]);
@@ -114 +114 @@
         /*get enum type of Penta: */
         enum_type=PentaEnum();
-        
+
         /*marshall enum: */
         memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
-        
+
         /*marshall Penta data: */
         memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
@@ -149 +149 @@
         memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
         memcpy(marshalled_dataset,&stokesreconditioning,sizeof(stokesreconditioning));marshalled_dataset+=sizeof(stokesreconditioning);
-        
+
         *pmarshalled_dataset=marshalled_dataset;
         return;
 }
-                
+
 int   Penta::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(friction_type)+
-                sizeof(p)+
-                sizeof(q)+
-                sizeof(shelf)+
-                sizeof(onbed)+
-                sizeof(onsurface)+
-                sizeof(meanvel)+
-                sizeof(epsvel)+
-                sizeof(collapse)+
-                sizeof(melting)+
-                sizeof(accumulation)+
-                sizeof(geothermalflux)+
-                sizeof(artdiff)+
-                sizeof(thermal_steadystate) +
-                sizeof(viscosity_overshoot) +
-                sizeof(stokesreconditioning) +
-                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(friction_type)+
+          sizeof(p)+
+          sizeof(q)+
+          sizeof(shelf)+
+          sizeof(onbed)+
+          sizeof(onsurface)+
+          sizeof(meanvel)+
+          sizeof(epsvel)+
+          sizeof(collapse)+
+          sizeof(melting)+
+          sizeof(accumulation)+
+          sizeof(geothermalflux)+
+          sizeof(artdiff)+
+          sizeof(thermal_steadystate) +
+          sizeof(viscosity_overshoot) +
+          sizeof(stokesreconditioning) +
+          sizeof(int); //sizeof(int) for enum type
 }
 
@@ -262 +262 @@
 
         int i;
-        
+
         DataSet* loadsin=NULL;
         DataSet* nodesin=NULL;
@@ -274 +274 @@
         /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */
         ResolvePointers((Object**)nodes,node_ids,node_offsets,6,nodesin);
-        
+
         /*Same for materials: */
         ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin);
@@ -293 +293 @@
         }
         else if (analysis_type==DiagnosticAnalysisEnum()){
-        
+
                 if (sub_analysis_type==HorizAnalysisEnum()){
-                
+
                         CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
                 }
                 else if (sub_analysis_type==VertAnalysisEnum()){
-                
+
                         CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type);
                 }
                 else if (sub_analysis_type==StokesAnalysisEnum()){
-                
+
                         CreateKMatrixDiagnosticStokes( Kgg,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()){
-                
+
                 CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==ThermalAnalysisEnum()){
-                
+
                 CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==MeltingAnalysisEnum()){
-                
+
                 CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type);
         }
@@ -342 +342 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-        
-        
+
+
         /* 3d gaussian points: */
         int     num_gauss,ig;
@@ -359 +359 @@
         int     num_area_gauss;
         double  gauss_weight;
-        
+
         /* 2d gaussian point: */
         int     num_gauss2d;
@@ -391 +391 @@
         double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
         double Jdet;
-        
+
         /*slope: */
         double  slope[2]={0.0,0.0};
@@ -454 +454 @@
                 }
 
-                #ifdef _DEBUGELEMENTS_
+#ifdef _DEBUGELEMENTS_
                 if(my_rank==RANK && id==ELID){ 
                         printf("El id %i Rank %i PentaElement input list before gaussian loop: \n",ELID,RANK); 
@@ -471 +471 @@
                         printf("   temperature_average [%g %g %g %g %g %g]\n",temperature_average_list[0],temperature_average_list[1],temperature_average_list[2],temperature_average_list[3],temperature_average_list[4],temperature_average_list[5]);
                 }
-                #endif
+#endif
 
                 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -484 +484 @@
                 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
 
-                #ifdef _DEBUGGAUSS_
+#ifdef _DEBUGGAUSS_
                 if(my_rank==RANK && id==ELID){ 
                         printf("El id %i Rank %i PentaElement gauss points\n",ELID,RANK); 
@@ -494 +494 @@
                         }       
                 }
-                #endif
+#endif
                 /* Start  looping on the number of gaussian points: */
                 for (ig1=0; ig1<num_area_gauss; ig1++){
                         for (ig2=0; ig2<num_vert_gauss; ig2++){
-                        
+
                                 /*Pick up the gaussian point: */
                                 gauss_weight1=*(area_gauss_weights+ig1);
                                 gauss_weight2=*(vert_gauss_weights+ig2);
                                 gauss_weight=gauss_weight1*gauss_weight2;
-                                
-                                
+
+
                                 gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
                                 gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
@@ -514 +514 @@
                                 GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
                                 GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
-                        
+
                                 /*Get viscosity: */
                                 matice->GetViscosity3d(&viscosity, &epsilon[0]);
@@ -525 +525 @@
                                 /* Get Jacobian determinant: */
                                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
-        
+
                                 /*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: */
-                                
+
                                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
                                 D_scalar=newviscosity*gauss_weight*Jdet;
@@ -534 +534 @@
                                         D[i][i]=D_scalar;
                                 }
-                
+
                                 /*  Do the triple product tB*D*Bprime: */
                                 TripleMultiply( &B[0][0],5,numdof,1,
-                                                &D[0][0],5,5,0,
-                                                &Bprime[0][0],5,numdof,0,
-                                                &Ke_gg_gaussian[0][0],0);
+                                                        &D[0][0],5,5,0,
+                                                        &Bprime[0][0],5,numdof,0,
+                                                        &Ke_gg_gaussian[0][0],0);
 
                                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
@@ -549 +549 @@
                         } //for (ig2=0; ig2<num_vert_gauss; ig2++)
                 } //for (ig1=0; ig1<num_area_gauss; ig1++)
-                
+
 
                 /*Add Ke_gg to global matrix Kgg: */
                 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
-        
+
                 //Deal with 2d friction at the bedrock interface
                 if((onbed && !shelf)){
@@ -567 +567 @@
 
         } 
-                
-        cleanup_and_return: 
+
+cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -627 +627 @@
         /*recover pointers: */
         inputs=(ParameterInputs*)vinputs;
-        
+
 
         /*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness 
@@ -660 +660 @@
 
         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-        #ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
         for (i=0;i<num_area_gauss;i++){
                 _printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
@@ -667 +667 @@
                 _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
         }
-        #endif
-        
+#endif
+
         /* Start  looping on the number of gaussian points: */
         for (ig1=0; ig1<num_area_gauss; ig1++){
                 for (ig2=0; ig2<num_vert_gauss; ig2++){
-                
+
                         /*Pick up the gaussian point: */
                         gauss_weight1=*(area_gauss_weights+ig1);
                         gauss_weight2=*(vert_gauss_weights+ig2);
                         gauss_weight=gauss_weight1*gauss_weight2;
-                        
+
                         gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
                         gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
@@ -693 +693 @@
                         /*  Do the triple product tB*D*Bprime: */
                         TripleMultiply( &B[0][0],1,numgrids,1,
-                                        &DL_scalar,1,1,0,
-                                        &Bprime[0][0],1,numgrids,0,
-                                        &Ke_gg_gaussian[0][0],0);
+                                                &DL_scalar,1,1,0,
+                                                &Bprime[0][0],1,numgrids,0,
+                                                &Ke_gg_gaussian[0][0],0);
 
                         /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
@@ -708 +708 @@
         /*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);
@@ -740 +740 @@
         double         gravity,rho_ice,rho_water;
 
-        
+
         /*Collapsed formulation: */
         Tria*  tria=NULL;
@@ -746 +746 @@
         /*Grid data: */
         double        xyz_list[numgrids][3];
-        
+
         /*parameters: */
         double             xyz_list_tria[numgrids2d][3];
@@ -770 +770 @@
         double     DLStokes[14][14]={0.0};
         double     tLDStokes[numdof2d][14];
-        
+
         /* gaussian points: */
         int     num_area_gauss;
@@ -792 +792 @@
         double  alpha2_list[numgrids2d];
         double  alpha2_gauss;
-        
+
         ParameterInputs* inputs=NULL;
 
@@ -804 +804 @@
                 }
         }
-        
+
         /*recovre material parameters: */
         rho_water=matpar->GetRhoWater();
@@ -818 +818 @@
 
         /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
-           get tria gaussian points as well as segment gaussian points. For tria gaussian 
-           points, order of integration is 2, because we need to integrate the product tB*D*B' 
-           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
-           points, same deal, which yields 3 gaussian points.*/
+                get tria gaussian points as well as segment gaussian points. For tria gaussian 
+                points, order of integration is 2, because we need to integrate the product tB*D*B' 
+                which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
+                points, same deal, which yields 3 gaussian points.*/
 
         area_order=5;
@@ -845 +845 @@
                         /*Compute strain rate: */
                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
-                
+
                         /*Get viscosity: */
                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
@@ -883 +883 @@
                 /*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();
@@ -910 +910 @@
                         }
                 }
-                
+
                 GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
 
@@ -920 +920 @@
                         gauss_coord[2]=*(third_gauss_area_coord+igarea);
                         gauss_coord[3]=-1;
-                        
+
                         gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); 
                         gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
                         gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
-        
+
                         /*Get the Jacobian determinant */
                         tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
-                        
+
                         /*Get L matrix if viscous basal drag present: */
                         GetLStokes(&LStokes[0][0],  gauss_coord_tria);
                         GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord);
-                                
+
                         /*Compute strain rate: */
                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
-                
+
                         /*Get viscosity at last iteration: */
                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
-        
+
                         /*Get normal vecyor to the bed */
                         SurfaceNormal(&surface_normal[0],xyz_list_tria);
-                        
+
                         bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
                         bed_normal[1]=-surface_normal[1];
@@ -974 +974 @@
                 }
         } //if ( (onbed==1) && (shelf==0))
-        
+
         /*Reduce the matrix */
         ReduceMatrixStokes(&Ke_reduced[0][0], &Ke_temp[0][0]);
@@ -986 +986 @@
         /*Add Ke_gg to global matrix Kgg: */
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
-        
-        cleanup_and_return: 
+
+cleanup_and_return: 
 
         return;
@@ -998 +998 @@
         /*Just branch to the correct element stiffness matrix 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()){
 
@@ -1008 +1008 @@
                 }
                 else if (sub_analysis_type==VertAnalysisEnum()){
-                
+
                         CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type);
                 }
                 else if (sub_analysis_type==StokesAnalysisEnum()){
-                
+
                         CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
                 }
@@ -1018 +1018 @@
         }
         else if (analysis_type==SlopeComputeAnalysisEnum()){
-                
+
                 CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==ThermalAnalysisEnum()){
-                
+
                 CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==MeltingAnalysisEnum()){
-                
+
                 CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type);
         }
@@ -1056 +1056 @@
         inputs->Recover("melting",&melting[0],1,dofs,6,(void**)nodes);
         inputs->Recover("accumulation_param",&accumulation[0],1,dofs,6,(void**)nodes);
-        
+
         //Update material if necessary
         if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){
@@ -1063 +1063 @@
                 matice->SetB(B_average);
         }
-        
+
         if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){
                 B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0;
@@ -1088 +1088 @@
         }
 }
-                
+
 int Penta::GetOnBed(){
         return onbed;
@@ -1099 +1099 @@
 }
 void          Penta::GetBedList(double* bed_list){
-        
+
         int i;
         for(i=0;i<6;i++)bed_list[i]=b[i];
@@ -1115 +1115 @@
         int i;
         Tria* tria=NULL;
-        
+
         /*Bail out if this element if:
          * -> Non collapsed and not on the surface
@@ -1132 +1132 @@
         }
         else{
-                
+
                 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
                 tria->Du(du_g,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
@@ -1156 +1156 @@
 #define __FUNCT__ "Penta::GradjDrag"
 void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){
-        
-        
+
+
         Tria* tria=NULL;
-        
+
         /*Bail out if this element does not touch the bedrock: */
         if (!onbed){
@@ -1165 +1165 @@
         }
         else{
-                
+
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
                 tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type);
@@ -1178 +1178 @@
         throw ErrorException(__FUNCT__," not supported yet!");
 }
-        
+
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::Misfit"
 double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type,int sub_analysis_type){
-        
+
         double J;
         Tria* tria=NULL;
-        
+
         /*Bail out if this element if:
          * -> Non collapsed and not on the surface
@@ -1209 +1209 @@
         }
 }
-                
+
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::SpawnTria"
@@ -1224 +1224 @@
         double   tria_accumulation[3];
         double   tria_geothermalflux[3];
-        
+
         /*configuration: */
         int tria_node_ids[3];
@@ -1249 +1249 @@
         tria_melting[1]=melting[g1];
         tria_melting[2]=melting[g2];
-        
+
         tria_accumulation[0]=accumulation[g0];
         tria_accumulation[1]=accumulation[g1];
@@ -1286 +1286 @@
         int doflist_per_node[MAXDOFSPERNODE];
         int numberofdofspernode;
-        
+
         for(i=0;i<6;i++){
                 nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
@@ -1320 +1320 @@
         GetB(&B[0][0], xyz_list, gauss_l1l2l3l4);
 
-        #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]);
@@ -1326 +1326 @@
         _printf_("              [ %lf   %lf  ]\n",B[3][0],B[3][1]);
         _printf_("              [ %lf   %lf  ]\n",B[4][0],B[4][1]);
-        
+
         _printf_("B for grid2 : [ %lf   %lf  \n",B[0][2],B[0][3]);
         _printf_("              [ %lf   %lf  \n",B[1][2],B[1][3]);
@@ -1332 +1332 @@
         _printf_("              [ %lf   %lf  ]\n",B[3][2],B[3][3]);
         _printf_("              [ %lf   %lf  ]\n",B[4][2],B[4][3]);
-        
+
         _printf_("B for grid3 : [ %lf   %lf  \n", B[0][4],B[0][5]);
         _printf_("              [ %lf   %lf  \n", B[1][4],B[1][5]);
@@ -1338 +1338 @@
         _printf_("              [ %lf   %lf  ]\n",B[3][4],B[3][5]);
         _printf_("              [ %lf   %lf  ]\n",B[4][4],B[4][5]);
-        
+
         _printf_("B for grid4 : [ %lf   %lf  \n", B[0][6],B[0][7]);
         _printf_("              [ %lf   %lf  \n", B[1][6],B[1][7]);
@@ -1344 +1344 @@
         _printf_("              [ %lf   %lf  ]\n",B[3][6],B[3][7]);
         _printf_("              [ %lf   %lf  ]\n",B[4][6],B[4][7]);
-                                
+
         _printf_("B for grid5 : [ %lf   %lf  \n", B[0][8],B[0][9]);
         _printf_("              [ %lf   %lf  \n", B[1][8],B[1][9]);
@@ -1360 +1360 @@
                 _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],5,NDOF2*numgrids,0,
-                                      velocity,NDOF2*numgrids,1,0,
-                                                  epsilon,0);
+                                velocity,NDOF2*numgrids,1,0,
+                                epsilon,0);
 
 }
@@ -1385 +1385 @@
          * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
          */
-        
+
         int i;
         const int numgrids=6;
@@ -1396 +1396 @@
         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
 
-        #ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
         for (i=0;i<numgrids;i++){
                 _printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
         }
-        #endif
+#endif
 
         /*Build B: */
@@ -1415 +1415 @@
                 *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i]; 
                 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
-                
+
                 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
                 *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i]; 
@@ -1439 +1439 @@
          * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
          */
-        
+
         int i;
         const int NDOF3=3;
@@ -1450 +1450 @@
         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
 
-        #ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
         for (i=0;i<numgrids;i++){
                 _printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
         }
-        #endif
+#endif
 
         /*Build BPrime: */
@@ -1469 +1469 @@
                 *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
                 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
-                
+
                 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
                 *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
@@ -1482 +1482 @@
          * the determinant of it: */
         const int NDOF3=3;
-        
+
         double J[NDOF3][NDOF3];
 
@@ -1496 +1496 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasic" 
 void Penta::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4){
-        
+
         /*This routine returns the values of the nodal functions derivatives  (with respect to the basic coordinate system: */
 
-        
+
         int i;
         const int NDOF3=3;
@@ -1534 +1534 @@
         const int NDOF3=3;
         int i,j;
-        
+
         /*The Jacobian is constant over the element, discard the gaussian points. 
          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
@@ -1544 +1544 @@
         double y1,y2,y3,y4,y5,y6;
         double z1,z2,z3,z4,z5,z6;
-        
+
         double sqrt3=sqrt(3.0);
-        
+
         /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
         A1=gauss_l1l2l3l4[0];
@@ -1562 +1562 @@
         x5=*(xyz_list+3*4+0);
         x6=*(xyz_list+3*5+0);
-        
+
         y1=*(xyz_list+3*0+1);
         y2=*(xyz_list+3*1+1);
@@ -1590 +1590 @@
         *(J+NDOF3*2+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+1.0/4.0*(z1-z2-z4+z5)*xi+1.0/4.0*(-z1+z5-z2+z4);
 
-        #ifdef _DEBUG_
+#ifdef _DEBUG_
         for(i=0;i<3;i++){
                 for (j=0;j<3;j++){
@@ -1597 +1597 @@
                 printf("\n");
         }
-        #endif
+#endif
 }
 
@@ -1603 +1603 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParams" 
 void Penta::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4){
-        
+
         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
          * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
@@ -1610 +1610 @@
         double A1,A2,A3,z;
         double sqrt3=sqrt(3.0);
-        
+
         A1=gauss_l1l2l3l4[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
         A2=gauss_l1l2l3l4[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
@@ -1677 +1677 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-        
+
         /* parameters: */
         double  slope[NDOF2];
@@ -1752 +1752 @@
 
                 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-                #ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
                 for (i=0;i<num_area_gauss;i++){
                         _printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
@@ -1759 +1759 @@
                         _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
                 }
-                #endif
-        
+#endif
+
                 /* Start  looping on the number of gaussian points: */
                 for (ig1=0; ig1<num_area_gauss; ig1++){
                         for (ig2=0; ig2<num_vert_gauss; ig2++){
-                        
+
                                 /*Pick up the gaussian point: */
                                 gauss_weight1=*(area_gauss_weights+ig1);
                                 gauss_weight2=*(vert_gauss_weights+ig2);
                                 gauss_weight=gauss_weight1*gauss_weight2;
-                                
+
                                 gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
                                 gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
                                 gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
                                 gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
-                
+
                                 /*Compute thickness at gaussian point: */
                                 GetParameterValue(&thickness, &h[0],gauss_l1l2l3l4);
@@ -1783 +1783 @@
                                 /* Get Jacobian determinant: */
                                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
-                
-                                 /*Get nodal functions: */
+
+                                /*Get nodal functions: */
                                 GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
 
@@ -1808 +1808 @@
         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);
@@ -1822 +1822 @@
 #define __FUNCT__ "Penta::CreateKMatrix"
 void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){
-        
+
         const int numgrids=6;
         double l1l2l3l4l5l6[numgrids];
@@ -1834 +1834 @@
 #define __FUNCT__ "Penta::GetParameterDerivativeValue"
 void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){
-                                
+
         /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_l1l2l3l4:
          *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
@@ -1849 +1849 @@
         /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
-        
+
         *(p+0)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[0][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[0][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[0][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[0][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[0][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[0][5];
-;
+        ;
         *(p+1)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[1][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[1][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[1][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[1][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[1][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[1][5];
 
@@ -1861 +1861 @@
 #define __FUNCT__ "Penta::GetNodalFunctions"
 void Penta::GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4){
-        
+
         /*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -1867 +1867 @@
 
         l1l2l3l4l5l6[1]=gauss_l1l2l3l4[1]*(1-gauss_l1l2l3l4[3])/2.0;
-        
+
         l1l2l3l4l5l6[2]=gauss_l1l2l3l4[2]*(1-gauss_l1l2l3l4[3])/2.0;
 
@@ -1873 +1873 @@
 
         l1l2l3l4l5l6[4]=gauss_l1l2l3l4[1]*(1+gauss_l1l2l3l4[3])/2.0;
-        
+
         l1l2l3l4l5l6[5]=gauss_l1l2l3l4[2]*(1+gauss_l1l2l3l4[3])/2.0;
 
@@ -1888 +1888 @@
         int          nodedofs[2];
         int          numberofdofspernode;
-        
+
         Node* node=NULL;
         int   i;
         double velocity[2];
-                
+
 
         if((collapse==1) && (onbed==1)){
-                        
+
                 GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1902 +1902 @@
                  * inserting the same velocity value into ug, until we reach the surface: */
                 for(i=0;i<3;i++){
-        
+
                         node=nodes[i]; //base nodes
-        
+
                         /*get velocity for this base node: */
                         velocity[0]=ug_serial[doflist[numberofdofspernode*i+0]];
@@ -1938 +1938 @@
         int          nodedof;
         int          numberofdofspernode;
-        
+
         Node* node=NULL;
         int   i;
         double slope;
-                
+
 
         if(onbed==1){
-                        
+
                 GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1952 +1952 @@
                 /*For each node on the base of this penta,  we grab the slope. Once we know the slope, we follow the upper nodes, 
                  * inserting the same slope value into sg, until we reach the surface: */
-                
+
                 for(i=0;i<3;i++){
-        
+
                         node=nodes[i]; //base nodes
-        
+
                         /*get velocity for this base node: */
                         slope=sg_serial[doflist[i]];
@@ -1984 +1984 @@
 
         /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
-        where hi is the interpolation function for grid i.*/
+                where hi is the interpolation function for grid i.*/
 
         int i;
@@ -1994 +1994 @@
         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
 
-        #ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
         for (i=0;i<numgrids;i++){
                 _printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
         }
-        #endif
+#endif
 
         /*Build B: */
@@ -2004 +2004 @@
                 B[i]=dh1dh2dh3dh4dh5dh6_basic[2][i];  
         }
-        
+
 }
 
@@ -2015 +2015 @@
 
         int i;
-                                
+
         GetNodalFunctions(B, gauss_l1l2l3l4);
 
@@ -2034 +2034 @@
         int          doflist[numdof];
         int          numberofdofspernode;
-        
+
         /* gaussian points: */
         int     num_gauss,ig;
@@ -2094 +2094 @@
         /* recover input parameters: */
         if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!");
-            inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
+        inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
 
         /*Get gaussian points and weights :*/
@@ -2101 +2101 @@
 
         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-        #ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
         for (i=0;i<num_area_gauss;i++){
                 _printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
@@ -2108 +2108 @@
                 _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
         }
-        #endif
+#endif
 
         /* Start  looping on the number of gaussian points: */
         for (ig1=0; ig1<num_area_gauss; ig1++){
                 for (ig2=0; ig2<num_vert_gauss; ig2++){
-                
+
                         /*Pick up the gaussian point: */
                         gauss_weight1=*(area_gauss_weights+ig1);
                         gauss_weight2=*(vert_gauss_weights+ig2);
                         gauss_weight=gauss_weight1*gauss_weight2;
-                        
+
                         gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
                         gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
                         gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
                         gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
-        
+
                         /*Get velocity derivative, with respect to x and y: */
                         GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
@@ -2129 +2129 @@
                         dudx=du[0];
                         dvdy=dv[1];
-                        
+
 
                         /* Get Jacobian determinant: */
                         GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
-                        #ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
                         _printf_("Element id %i Jacobian determinant: %lf\n",PentaElementGetID(this),Jdet);
-                        #endif
-                
-                         /*Get nodal functions: */
+#endif
+
+                        /*Get nodal functions: */
                         GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
 
@@ -2154 +2154 @@
         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);
@@ -2175 +2175 @@
         double rho_ice,g;
         double       xyz_list[numgrids][3];
-                
+
         /*Get node data: */
         GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
-        
+
         /*pressure is lithostatic: */
         //md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab
@@ -2191 +2191 @@
                 pressure[i]=rho_ice*g*(s[i]-xyz_list[i][2]);
         }
-        
+
         /*plug local pressure values into global pressure vector: */
         VecSetValues(pg,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
@@ -2205 +2205 @@
         /*Collapsed formulation: */
         Tria*  tria=NULL;
-        
+
         /*Is this element on the bed? :*/
         if(!onbed)return;
@@ -2221 +2221 @@
 
 void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
-        
+
         /*Collapsed formulation: */
         Tria*  tria=NULL;
-        
+
         /*Is this element on the bed? :*/
         if(!onbed)return;
@@ -2238 +2238 @@
 #define __FUNCT__ "ReduceMatrixStokes" 
 void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
-                
+
         int i,j;
 
@@ -2292 +2292 @@
         double det;
         int calculationdof=3;
-        
+
         /*Take the matrix components: */
         a=*(Ke+calculationdof*0+0);
@@ -2305 +2305 @@
 
         det=a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g);
-        
+
         *(Ke_invert+calculationdof*0+0)=(e*i-f*h)/det;
         *(Ke_invert+calculationdof*0+1)=(c*h-b*i)/det;
@@ -2392 +2392 @@
          * For grid i, Bi can be expressed in the basic coordinate system
          * by:                          Bi_basic=[ dh/dx          0             0       0  ]
-                *                                       [   0           dh/dy           0       0  ]
-                *                                       [   0             0           dh/dy     0  ]
-                *                                       [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
-                *                                       [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
-                *                                       [   0          1/2*dh/dz    1/2*dh/dy   0  ]
-                *                                       [   0             0             0       h  ]
-                *                                       [ dh/dx         dh/dy         dh/dz     0  ]
-                *       where h is the interpolation function for grid i.
-                *       Same thing for Bb except the last column that does not exist.
-                */
-        
+         *                                      [   0           dh/dy           0       0  ]
+         *                                      [   0             0           dh/dy     0  ]
+         *                                      [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
+         *                                      [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
+         *                                      [   0          1/2*dh/dz    1/2*dh/dy   0  ]
+         *                                      [   0             0             0       h  ]
+         *                                      [ dh/dx         dh/dy         dh/dz     0  ]
+         *      where h is the interpolation function for grid i.
+         *      Same thing for Bb except the last column that does not exist.
+         */
+
         int i;
         int calculationdof=3;
@@ -2416 +2416 @@
 
         GetNodalFunctions(l1l6, gauss_coord);
-        
-        #ifdef _DEBUG_ 
+
+#ifdef _DEBUG_ 
         for (i=0;i<7;i++){
                 printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i],dh1dh7_basic[2][i]);
         }
 
-        #endif
+#endif
 
         /*Build B: */
@@ -2470 +2470 @@
 
         /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 
-        *       For grid i, Bi' can be expressed in the basic coordinate system
-        *       by: 
-        *                               Bi_basic'=[ dh/dx   0          0       0]
-        *                                        [   0      dh/dy      0       0]
-        *                                        [   0      0         dh/dz    0]
-        *                                        [  dh/dy   dh/dx      0       0]
-        *                                        [  dh/dz   0        dh/dx     0]
-        *                                        [   0      dh/dz    dh/dy     0]
-        *                                        [  dh/dx   dh/dy    dh/dz     0]
-        *                                        [   0      0          0       h]
-        *       where h is the interpolation function for grid i.
-        *
-        *       Same thing for the bubble fonction except that there is no fourth column
-        */
+         *      For grid i, Bi' can be expressed in the basic coordinate system
+         *      by: 
+         *                              Bi_basic'=[ dh/dx   0          0       0]
+         *                                       [   0      dh/dy      0       0]
+         *                                       [   0      0         dh/dz    0]
+         *                                       [  dh/dy   dh/dx      0       0]
+         *                                       [  dh/dz   0        dh/dx     0]
+         *                                       [   0      dh/dz    dh/dy     0]
+         *                                       [  dh/dx   dh/dy    dh/dz     0]
+         *                                       [   0      0          0       h]
+         *      where h is the interpolation function for grid i.
+         *
+         *      Same thing for the bubble fonction except that there is no fourth column
+         */
 
         int i;
@@ -2498 +2498 @@
 
         GetNodalFunctions(l1l6, gauss_coord);
-        
-        #ifdef _DEBUG_ 
+
+#ifdef _DEBUG_ 
         for (i=0;i<6;i++){
                 printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
         }
 
-        #endif
+#endif
 
         /*B_primeuild B_prime: */
@@ -2551 +2551 @@
 
         /*
-        * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-        * For grid i, Li can be expressed in the basic coordinate system
-        * by: 
-        *       Li_basic=[ h    0    0   0]
-        *                [ 0    h    0   0]
-        *                [ 0    0    h   0]
-        *                [ 0    0    h   0]
-        *                [ h    0    0   0]
-        *                [ 0    h    0   0]
-        *                [ h    0    0   0]
-        *                [ 0    h    0   0]
-        *                [ 0    0    h   0]
-        *                [ 0    0    h   0]
-        *                [ 0    0    h   0]
-        *                [ h    0    0   0]
-        *                [ 0    h    0   0]
-        *                [ 0    0    h   0]
-        * where h is the interpolation function for grid i.
-        */
+         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+         * For grid i, Li can be expressed in the basic coordinate system
+         * by: 
+         *       Li_basic=[ h    0    0   0]
+         *               [ 0    h    0   0]
+         *               [ 0    0    h   0]
+         *               [ 0    0    h   0]
+         *               [ h    0    0   0]
+         *               [ 0    h    0   0]
+         *               [ h    0    0   0]
+         *               [ 0    h    0   0]
+         *               [ 0    0    h   0]
+         *               [ 0    0    h   0]
+         *               [ 0    0    h   0]
+         *               [ h    0    0   0]
+         *               [ 0    h    0   0]
+         *               [ 0    0    h   0]
+         * where h is the interpolation function for grid i.
+         */
 
         int i;
@@ -2582 +2582 @@
         l1l2l3[1]=gauss_coord_tria[1];
         l1l2l3[2]=gauss_coord_tria[2];
-        
-
-        #ifdef _DELUG_ 
+
+
+#ifdef _DELUG_ 
         for (i=0;i<3;i++){
                 printf("Node %i  h=%lf \n",i,l1l2l3[i]);
         }
-        #endif
+#endif
 
         /*Build LStokes: */
@@ -2648 +2648 @@
                 *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
                 *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
-        
+
         }
 }
@@ -2658 +2658 @@
 
         /*
-        * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
-        * For grid i, Lpi can be expressed in the basic coordinate system
-        * by: 
-        *       Lpi_basic=[ h    0    0   0]
-        *                [ 0    h    0   0]
-        *                [ h    0    0   0]
-        *                [ 0    h    0   0]
-        *                [ 0    0    h   0]
-        *                [ 0    0    h   0]
-        *                [ 0    0  dh/dz 0]
-        *                [ 0    0  dh/dz 0]
-        *                [ 0    0  dh/dz 0]
-        *                [dh/dz 0  dh/dx 0]
-        *                [ 0 dh/dz dh/dy 0]
-        *                [ 0    0    0   h]
-        *                [ 0    0    0   h]
-        *                [ 0    0    0   h]
-        * where h is the interpolation function for grid i.
-        */
+         * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
+         * For grid i, Lpi can be expressed in the basic coordinate system
+         * by: 
+         *       Lpi_basic=[ h    0    0   0]
+         *               [ 0    h    0   0]
+         *               [ h    0    0   0]
+         *               [ 0    h    0   0]
+         *               [ 0    0    h   0]
+         *               [ 0    0    h   0]
+         *               [ 0    0  dh/dz 0]
+         *               [ 0    0  dh/dz 0]
+         *               [ 0    0  dh/dz 0]
+         *               [dh/dz 0  dh/dx 0]
+         *               [ 0 dh/dz dh/dy 0]
+         *               [ 0    0    0   h]
+         *               [ 0    0    0   h]
+         *               [ 0    0    0   h]
+         * where h is the interpolation function for grid i.
+         */
 
         int i;
@@ -2690 +2690 @@
         l1l2l3[1]=gauss_coord_tria[1];
         l1l2l3[2]=gauss_coord_tria[2];
-        
+
         GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
 
-        #ifdef _DELUG_ 
+#ifdef _DELUG_ 
         for (i=0;i<3;i++){
                 printf("Node %i  h=%lf \n",i,l1l2l3[i]);
         }
-        #endif
+#endif
 
         /*Build LprimeStokes: */
@@ -2757 +2757 @@
                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
-        
-        }
-        
+
+        }
+
 }
 
@@ -2765 +2765 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasicStokes"
 void Penta::GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord){
-        
+
         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
          * basic coordinate system: */
@@ -2801 +2801 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParamsStokes"
 void Penta::GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord){
-        
+
         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
          * natural coordinate system) at the gaussian point. */
@@ -2893 +2893 @@
         int     area_order=5;
         int       num_vert_gauss=5;
-        
+
         double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         double  viscosity;
@@ -2916 +2916 @@
         double     tBD[27][8];
         double     P_terms[numdof];
-        
+
         ParameterInputs* inputs=NULL;
         Tria*            tria=NULL;
@@ -2941 +2941 @@
 
         /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
-           get tria gaussian points as well as segment gaussian points. For tria gaussian 
-           points, order of integration is 2, because we need to integrate the product tB*D*B' 
-           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
-           points, same deal, which yields 3 gaussian points.*/
+                get tria gaussian points as well as segment gaussian points. For tria gaussian 
+                points, order of integration is 2, because we need to integrate the product tB*D*B' 
+                which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
+                points, same deal, which yields 3 gaussian points.*/
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -2964 +2964 @@
                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
-                
+
                         /* Get Jacobian determinant: */
                         GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
@@ -2984 +2984 @@
                         GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); 
                         GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); 
-                        
+
                         /*Get bubble part of Bprime */
                         for(i=0;i<8;i++){
@@ -3033 +3033 @@
                         gauss_coord[2]=*(third_gauss_area_coord+igarea);
                         gauss_coord[3]=-1;
-                        
+
                         gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); 
                         gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
@@ -3040 +3040 @@
                         /*Get the Jacobian determinant */
                         tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
-                        
+
                         /* Get bed at gaussian point */
                         GetParameterValue(&bed,&b[0],gauss_coord);
@@ -3046 +3046 @@
                         /*Get L matrix : */
                         tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1);
-                                
+
                         /*Get water_pressure at gaussian point */
                         water_pressure=gravity*rho_water*bed;
@@ -3052 +3052 @@
                         /*Get normal vecyor to the bed */
                         SurfaceNormal(&surface_normal[0],xyz_list_tria);
-                        
+
                         bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
                         bed_normal[1]=-surface_normal[1];
@@ -3064 +3064 @@
                 }
         } //if ( (onbed==1) && (shelf==1))
-        
+
         /*Reduce the matrix */
         ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]);
@@ -3080 +3080 @@
 #define __FUNCT__ "Penta::ReduceVectorStokes" 
 void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){
-                
+
         int i,j;
 
@@ -3123 +3123 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsStokes" 
 void Penta::GetNodalFunctionsStokes(double* l1l7, double* gauss_coord){
-        
+
         /*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -3161 +3161 @@
         const int    NDOF1=1;
         const int    numdof=NDOF1*numgrids;
-        double       xyz_list[numgrids][3];
-        int          doflist[numdof];
-        int          numberofdofspernode;
+        double xyz_list[numgrids][3];
+        int    doflist[numdof];
+        int    numberofdofspernode;
 
         /* gaussian points: */
@@ -3177 +3177 @@
 
         int area_order=5;
-        int     num_vert_gauss=5;
-        
-        int            dofs[3]={0,1,2};
-        double         dt;
-        int            dt_exists;
-        
-        double         vxvyvz_list[numgrids][3];
-        double         vx_list[numgrids];
-        int            vx_list_exists;
-        double         u;
-        double         vy_list[numgrids];
-        int            vy_list_exists;
-        double         v;
-        double         vz_list[numgrids];
-        int            vz_list_exists;
-        double         w;
-
-        /* element parameters: */
-        int         onbed;
-        int         shelf;
-        int     steady_state;
+        int num_vert_gauss=5;
+
+        int     dofs[3]={0,1,2};
+        double  dt;
+        int     dt_exists;
+
+        double  vxvyvz_list[numgrids][3];
+        double  vx_list[numgrids];
+        int     vx_list_exists;
+        double  u;
+        double  vy_list[numgrids];
+        int     vy_list_exists;
+        double  v;
+        double  vz_list[numgrids];
+        int     vz_list_exists;
+        double  w;
 
         /*matrices: */
-        double     K_terms[numdof][numdof]={0.0};
-        double     Ke_gaussian_conduct[numdof][numdof];
-        double     Ke_gaussian_advec[numdof][numdof];
-        double     Ke_gaussian_transient[numdof][numdof];
-        double     B[3][numdof];
-        double     Bprime[3][numdof];
-        double     B_conduct[3][numdof];
-        double     B_advec[3][numdof];
-        double     Bprime_advec[3][numdof];
-        double     L[numdof];
-        double     D_scalar;
-        double     D[3][3];
-        double     l1l2l3[3];
-        double     tl1l2l3D[3];
-        double     tBD[3][numdof];
-        double     tBD_conduct[3][numdof];
-        double     tBD_advec[3][numdof];
-        double     tLD[numdof];
-
-        double     Jdet;
+        double  K_terms[numdof][numdof]={0.0};
+        double  Ke_gaussian_conduct[numdof][numdof];
+        double  Ke_gaussian_advec[numdof][numdof];
+        double  Ke_gaussian_transient[numdof][numdof];
+        double  B[3][numdof];
+        double  Bprime[3][numdof];
+        double  B_conduct[3][numdof];
+        double  B_advec[3][numdof];
+        double  Bprime_advec[3][numdof];
+        double  L[numdof];
+        double  D_scalar;
+        double  D[3][3];
+        double  l1l2l3[3];
+        double  tl1l2l3D[3];
+        double  tBD[3][numdof];
+        double  tBD_conduct[3][numdof];
+        double  tBD_advec[3][numdof];
+        double  tLD[numdof];
+
+        double  Jdet;
 
         /*Material properties: */
-        double         gravity,rho_ice,rho_water;
-        double         heatcapacity,thermalconductivity;
-        double         mixed_layer_capacity,thermal_exchange_velocity;
+        double  gravity,rho_ice,rho_water;
+        double  heatcapacity,thermalconductivity;
+        double  mixed_layer_capacity,thermal_exchange_velocity;
 
         /*Collapsed formulation: */
@@ -3237 +3232 @@
         GetDofList(&doflist[0],&numberofdofspernode);
         
-        // /*recovre material parameters: */
+        /*recovre material parameters: */
         rho_water=matpar->GetRhoWater();
         rho_ice=matpar->GetRhoIce();
@@ -3243 +3238 @@
         heatcapacity=matpar->GetHeatCapacity();
         thermalconductivity=matpar->GetThermalConductivity();
-
                 
         /*recover extra inputs from users, dt and velocity: */
@@ -3300 +3294 @@
                         D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
 
-
                         /*  Do the triple product B'*D*B: */
                         MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0);
@@ -3310 +3303 @@
                         GetB_advec(&B_advec[0][0],&xyz_list[0][0],gauss_coord); 
                         GetBprime_advec(&Bprime_advec[0][0],&xyz_list[0][0],gauss_coord); 
-
 
                         //Build the D matrix
@@ -3330 +3322 @@
                         MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0);
                         MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0);
-
 
                         /*Transient: */
@@ -3359 +3350 @@
                 }
         }
-        
-        //Ice/ocean heat exchange flux on ice shelf base 
-
-        if(onbed && shelf){
-
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type);
-                delete tria;
-                return;
-        }
 
         /*Add Ke_gg to global matrix Kgg: */
@@ -3381 +3362 @@
         xfree((void**)&vert_gauss_coord);
 
+        //Ice/ocean heat exchange flux on ice shelf base 
+        if(onbed && shelf){
+
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->CreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type);
+                delete tria;
+        }
 }
 
@@ -3397 +3385 @@
          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
          */
-        
+
         int i;
         int calculationdof=3;
@@ -3431 +3419 @@
          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
          */
-        
+
         int i;
         int calculationdof=3;
@@ -3466 +3454 @@
          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
          */
-        
+
         int i;
         int calculationdof=3;
@@ -3489 +3477 @@
 #define __FUNCT__ "Penta::CreateKMatrixMelting"
 void  Penta::CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
-        
+
         Tria* tria=NULL;
         if (!onbed){
@@ -3495 +3483 @@
         }
         else{
-                
+
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
                 tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type);
@@ -3520 +3508 @@
         /*Grid data: */
         double        xyz_list[numgrids][3];
-        
+
         /* gaussian points: */
         int     num_area_gauss,igarea,igvert;
@@ -3533 +3521 @@
         int area_order=2;
         int     num_vert_gauss=3;
-        
+
         double         dt;
         double         vx_list[numgrids];
@@ -3549 +3537 @@
         /* element parameters: */
         int     friction_type;
-        
+
         int            dofs[3]={0,1,2};
         int            dofs1[1]={0};
@@ -3584 +3572 @@
         GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-        
+
         // /*recovre material parameters: */
         rho_water=matpar->GetRhoWater();
@@ -3604 +3592 @@
                 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
         }
-        
+
         for(i=0;i<numgrids;i++){
                 vx_list[i]=vxvyvz_list[i][0];
@@ -3612 +3600 @@
 
         /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
-           get tria gaussian points as well as segment gaussian points. For tria gaussian 
-           points, order of integration is 2, because we need to integrate the product tB*D*B' 
-           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
-           points, same deal, which yields 3 gaussian points.: */
-        
+                get tria gaussian points as well as segment gaussian points. For tria gaussian 
+                points, order of integration is 2, because we need to integrate the product tB*D*B' 
+                which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
+                points, same deal, which yields 3 gaussian points.: */
+
         /*Get gaussian points: */
         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
@@ -3631 +3619 @@
                         gauss_coord[2]=*(third_gauss_area_coord+igarea);
                         gauss_coord[3]=*(vert_gauss_coord+igvert);
-        
+
                         /*Compute strain rate and viscosity: */
                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
@@ -3673 +3661 @@
         /* Ice/ocean heat exchange flux on ice shelf base */
         if(onbed && shelf){
-
 
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.

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

@@ -16 +16 @@
 #include "../DataSet/DataSet.h"
 #include "../include/typedefs.h"
-
 
 /*For debugging purposes: */
@@ -51 +50 @@
                 melting[i]=tria_melting[i];
                 accumulation[i]=tria_accumulation[i];
-                geothermalflux[i]=tria_geothermalflux[i]; }
+                geothermalflux[i]=tria_geothermalflux[i];
+        }
         matice=NULL;
         matice_offset=UNDEF;
@@ -242 +242 @@
 }
 
-
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::Configure"
@@ -294 +293 @@
                 throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
         }
-
-}
-
+}
 
 #undef __FUNCT__ 
@@ -302 +299 @@
 
 void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
-
 
         /* local declarations */