Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5741)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5742)
@@ -867,11 +867,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_coord[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	double  surface_list[5];
@@ -911,7 +906,4 @@
 	nx[0]=normal1[0]; nx[1]=normal2[0]; nx[2]=normal3[0]; nx[3]=normal4[0];
 	ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1];
-
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 	
 	//Recover the surface of the four nodes
@@ -976,20 +968,15 @@
 	xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45);
 
-	
-
-	//Start  looping on the triangle's gaussian points: 
-	for(ig=0;ig<num_gauss;ig++){
-
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_coord[0]=*(first_gauss_area_coord+ig);
-		gauss_coord[1]=*(second_gauss_area_coord+ig);
-		gauss_coord[2]=*(third_gauss_area_coord+ig);
-        
-		tria->GetNodalFunctions(l1l2l3_tria,gauss_coord);
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for(ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		tria->GetNodalFunctions(l1l2l3_tria,gauss);
 
 		//Get the coordinates of gauss point for each triangle in penta/quad
-		r_tria=gauss_coord[1]-gauss_coord[0];
-		s_tria=-3/sqrt((double)3)*(gauss_coord[0]+gauss_coord[1]-2/3);
+		r_tria=gauss->coord2-gauss->coord1;
+		s_tria=-3/sqrt((double)3)*(gauss->coord1+gauss->coord2-2/3);
 
 		//Coordinates of gauss points in the reference triangle
@@ -1070,26 +1057,20 @@
 
 
-			pe_g[0]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i];
-			pe_g[1]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i];
-			pe_g[2]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i];
-			pe_g[3]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i];
-			pe_g[4]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i];
-			pe_g[5]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i];
-			pe_g[6]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i];
-			pe_g[7]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i];
+			pe_g[0]+= J[i]*gauss->weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i];
+			pe_g[1]+= J[i]*gauss->weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i];
+			pe_g[2]+= J[i]*gauss->weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i];
+			pe_g[3]+= J[i]*gauss->weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i];
+			pe_g[4]+= J[i]*gauss->weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i];
+			pe_g[5]+= J[i]*gauss->weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i];
+			pe_g[6]+= J[i]*gauss->weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i];
+			pe_g[7]+= J[i]*gauss->weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i];
 
 			
 
 		} //for(i=0;i<4;i++)
-	} //for(ig=0;ig<num_gauss;ig++)
-	
-	/*Free ressources: */
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
-
-	/*Delete fake tria: */
+	}
+	
 	delete tria;
+	delete gauss;
 }
 /*}}}*/
@@ -1129,11 +1110,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_coord[3];
+	int     ig;
+	GaussTria* gauss=NULL;
 
 	double  surface_list[5];
@@ -1172,7 +1148,4 @@
 	ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1];
 	nz[0]=normal1[2]; nz[1]=normal2[2]; nz[2]=normal3[2]; nz[3]=normal4[2];
-
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 	
 	//Recover the surface of the four nodes
@@ -1237,20 +1210,15 @@
 	xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45);
 
-	
-
-	//Start  looping on the triangle's gaussian points: 
-	for(ig=0;ig<num_gauss;ig++){
-
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_coord[0]=*(first_gauss_area_coord+ig);
-		gauss_coord[1]=*(second_gauss_area_coord+ig);
-		gauss_coord[2]=*(third_gauss_area_coord+ig);
-        
-		tria->GetNodalFunctions(l1l2l3_tria,gauss_coord);
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for(ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		tria->GetNodalFunctions(l1l2l3_tria,gauss);
 
 		//Get the coordinates of gauss point for each triangle in penta/quad
-		r_tria=gauss_coord[1]-gauss_coord[0];
-		s_tria=-3/sqrt((double)3)*(gauss_coord[0]+gauss_coord[1]-2/3);
+		r_tria=gauss->coord2-gauss->coord1;
+		s_tria=-3/sqrt((double)3)*(gauss->coord1+gauss->coord2-2/3);
 
 		//Coordinates of gauss points in the reference triangle
@@ -1324,32 +1292,26 @@
 			pressure_tria = water_pressure_tria + air_pressure_tria;
 
-			pe_g[0]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*nx[i];
-			pe_g[1]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*ny[i];
-			pe_g[2]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*nz[i];
+			pe_g[0]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*nx[i];
+			pe_g[1]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*ny[i];
+			pe_g[2]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*nz[i];
 			pe_g[3]+= 0;
-			pe_g[4]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*nx[i];
-			pe_g[5]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*ny[i];
-			pe_g[6]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*nz[i];
+			pe_g[4]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*nx[i];
+			pe_g[5]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*ny[i];
+			pe_g[6]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*nz[i];
 			pe_g[7]+= 0;
-			pe_g[8]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*nx[i];
-			pe_g[9]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*ny[i];
-			pe_g[10]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*nz[i];
+			pe_g[8]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*nx[i];
+			pe_g[9]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*ny[i];
+			pe_g[10]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*nz[i];
 			pe_g[11]+= 0;
-			pe_g[12]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*nx[i];
-			pe_g[13]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*ny[i];
-			pe_g[14]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*nz[i];
+			pe_g[12]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*nx[i];
+			pe_g[13]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*ny[i];
+			pe_g[14]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*nz[i];
 			pe_g[15]+= 0;
 
 		} //for(i=0;i<4;i++)
 	} //for(ig=0;ig<num_gauss;ig++)
-	
-	/*Free ressources: */
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
-
-	/*Delete fake tria: */
+
 	delete tria;
+	delete gauss;
 }
 /*}}}*/