Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5623)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5624)
@@ -563,5 +563,5 @@
 
 	/* 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);
+	GaussLegendreTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2);
 
 	/*Retrieve all inputs we will be needing: */
@@ -2499,5 +2499,5 @@
 	  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 
+	  which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian 
 	  points, same deal, which yields 3 gaussian points.*/
 
@@ -2505,5 +2505,5 @@
 	num_vert_gauss=5;
 
-	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);
+	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);
 
 	/*Retrieve all inputs we will be needing: */
@@ -2830,5 +2830,5 @@
 	  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 
+	  which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian 
 	  points, same deal, which yields 3 gaussian points.*/
 
@@ -2836,5 +2836,5 @@
 	num_vert_gauss=5;
 
-	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);
+	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);
 
 	/*Retrieve all inputs we will be needing: */
@@ -3021,5 +3021,5 @@
 		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 
+		which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian 
 		points, same deal, which yields 3 gaussian points.*/
 
@@ -3028,5 +3028,5 @@
 
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	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);
+	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);
 
 	/*Retrieve all inputs we will be needing: */
@@ -3097,5 +3097,5 @@
 
 		xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
-		GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
+		GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
 
 		/* Start looping on the number of gauss 2d (nodes on the bedrock) */
@@ -3249,5 +3249,5 @@
 	  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 
+	  which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian 
 	  points, same deal, which yields 3 gaussian points.*/
 
@@ -3255,5 +3255,5 @@
 	num_vert_gauss=2;
 
-	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);
+	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);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3502,5 +3502,5 @@
 		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 
+		which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian 
 		points, same deal, which yields 3 gaussian points.: */
 
@@ -3509,5 +3509,5 @@
 	num_vert_gauss=2;
 
-	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);
+	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);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3861,5 +3861,5 @@
 	//Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
 	num_gauss=3;
-	GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
+	gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
 
 	/*Loop on the three segments*/
@@ -4034,5 +4034,5 @@
 	num_vert_gauss=3;
 
-	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);
+	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);
 
 	/*Retrieve all inputs we will be needing: */
@@ -4197,9 +4197,9 @@
 		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 
+		which is a polynomial of degree 3 (see GaussLegendreTria 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): */
-	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);
+	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);
 
 	/*Retrieve all inputs we will be needing: */
@@ -4285,5 +4285,5 @@
 
 		xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
-		GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
+		GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
 
 		/* Start looping on the number of gauss 2d (nodes on the bedrock) */
@@ -4437,5 +4437,5 @@
 	num_vert_gauss=2;
 
-	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);
+	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);
 
 	/*Retrieve all inputs we will be needing: */
@@ -4661,9 +4661,9 @@
 		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 
+		which is a polynomial of degree 3 (see GaussLegendreTria 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);
+	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);
 
 	/*Retrieve all inputs we will be needing: */
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5623)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5624)
@@ -727,5 +727,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -1081,5 +1081,5 @@
 
 	/* 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, 4);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
 
 	/*Retrieve all inputs we will be needing: */
@@ -1254,5 +1254,5 @@
 
 	/* 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, 4);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
 
 	/*Retrieve all inputs we will be needing: */
@@ -2003,5 +2003,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2128,5 +2128,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2260,5 +2260,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2390,5 +2390,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2522,5 +2522,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2653,5 +2653,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3105,5 +3105,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3232,5 +3232,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needed*/
@@ -3396,5 +3396,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3537,5 +3537,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needing: */
@@ -3695,5 +3695,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3838,5 +3838,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3979,5 +3979,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -4118,5 +4118,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Build normal vector to the surface:*/
@@ -4232,5 +4232,5 @@
 
 	/* Get gaussian points and weights: */
-	GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+	GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start looping on the number of gauss  (nodes on the bedrock) */
@@ -4358,5 +4358,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -4506,5 +4506,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needed*/
@@ -4614,5 +4614,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -4705,5 +4705,5 @@
 	heatcapacity=matpar->GetHeatCapacity();
 
-	GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+	GaussLegendreTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start looping on the number of gauss (nodes on the bedrock) */
@@ -4790,5 +4790,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*retrieve inputs :*/
@@ -4874,5 +4874,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needing: */
@@ -4956,5 +4956,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needing: */
@@ -5047,5 +5047,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needing: */
@@ -5166,5 +5166,5 @@
 
 	/* Get gaussian points and weights: */
-	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*/
+	GaussLegendreTria( &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*/
 
 	/* Start  looping on the number of gaussian points: */
@@ -5275,5 +5275,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -5511,5 +5511,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -5750,5 +5750,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start  looping on the number of gaussian points: */
@@ -5908,5 +5908,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needing: */
@@ -5996,5 +5996,5 @@
 
 	/* 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);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needing: */
@@ -6082,5 +6082,5 @@
 
 	/* Get gaussian points and weights: */
-	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*/
+	GaussLegendreTria( &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*/
 
 	/*Retrieve all inputs we will be needing: */
@@ -6193,5 +6193,5 @@
 	/* Ice/ocean heat exchange flux on ice shelf base */
 
-	GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+	GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/*Retrieve all inputs we will be needing: */
@@ -6325,5 +6325,5 @@
 	
 	/* Ice/ocean heat exchange flux on ice shelf base */
-	GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+	GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
 
 	/* Start looping on the number of gauss 2d (nodes on the bedrock) */
@@ -6844,5 +6844,5 @@
 
 	/* 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, 4);
+	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
 
 	/* Start  looping on the number of gaussian points: */
Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5623)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5624)
@@ -456,5 +456,5 @@
 	//Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
 	num_gauss=3;
-	GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
+	gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
 
 	for(ig=0;ig<num_gauss;ig++){
@@ -916,5 +916,5 @@
 
 	/* 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);
+	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
@@ -1177,5 +1177,5 @@
 
 	/* 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);
+	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
Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp	(revision 5623)
+++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp	(revision 5624)
@@ -455,5 +455,5 @@
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
 	num_gauss=2;
-	GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
+	gaussSegment(&gauss_coords, &gauss_weights,num_gauss);
 
 	/* Start  looping on the number of gaussian points: */
@@ -598,5 +598,5 @@
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
 	num_gauss=2;
-	GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
+	gaussSegment(&gauss_coords, &gauss_weights,num_gauss);
 
 	/* Start  looping on the number of gaussian points: */
@@ -745,5 +745,5 @@
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
 	num_gauss=2;
-	GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
+	gaussSegment(&gauss_coords, &gauss_weights,num_gauss);
 
 	/* Start  looping on the number of gaussian points: */
Index: /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp	(revision 5623)
+++ /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp	(revision 5624)
@@ -10,6 +10,6 @@
 
 /*General Gauss points*/
-/*FUNCTION GaussLegendre {{{1*/
-void GaussLegendre( double** pxgaus, double** pxwgt, int ngaus){
+/*FUNCTION GaussLegendreLinear {{{1*/
+void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus){
 	/* Gauss-Legendre quadrature points.
 
@@ -88,246 +88,6 @@
 	}
 }/*}}}1*/
-/*FUNCTION GaussLobatto{{{1*/
-void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {
-	/*Gauss-Lobatto quadrature points.
-
-	  The recurrence coefficients for Legendre polynomials on (-1,1)
-	  are defined (from the ORTHPOL subroutine RECUR with ipoly=1) as:
-
-	  alpha(i)=0.
-	  beta (i)=1./(4.-1./(i-1)^2))
-
-	  and then modified for the Gauss-Lobatto quadrature rule on (-1,1)
-	  (from the ORTHPOL subroutine LOB).
-
-	  For degree p, the required number of Gauss-Lobatto points is
-	  n>=(p+1)/2+1 (one more than Gauss-Legendre).*/
-
-	int i;
-	double *alpha,*beta;
-	double left=-1.,right= 1.;
-	double p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;
-
-	/*p= 1, npoint= 1 (Gauss-Legendre)*/
-	static double wgt1[]={2.000000000000000};
-	static double xi1[]={0.000000000000000};
-
-	/*p= 1, npoint= 2*/
-	static double wgt2[]={1.000000000000000, 1.000000000000000};
-	static double xi2[]={-1.000000000000000, 1.000000000000000};
-
-	/*p= 3, npoint= 3*/
-	static double wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
-	static double xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
-
-	/*p= 5, npoint= 4*/
-	static double wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};
-	static double xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};
-
-	/*p= 7, npoint= 5*/
-	static double wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};
-	static double xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};
-
-	static double* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };
-	static double* xip [MAX_LINE_GLOB_PTS]={xi1  ,xi2  ,xi3  ,xi4  ,xi5  };
-
-	static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof(double),
-		sizeof(wgt2 )/sizeof(double),
-		sizeof(wgt3 )/sizeof(double),
-		sizeof(wgt4 )/sizeof(double),
-		sizeof(wgt5 )/sizeof(double)};
-
-	//	_printf_("Gauss-Lobatto recurrence coefficients ngaus=%d\n",ngaus);
-	*pxgaus = (double *) xmalloc(ngaus*sizeof(double));
-	*pxwgt  = (double *) xmalloc(ngaus*sizeof(double));
-
-	/*  check to see if Gauss points need to be calculated  */
-	if (ngaus <= MAX_LINE_GLOB_PTS) {
-
-		/*  copy the points from the static arrays (noting that the pointers
-			 could be returned directly, but then the calling function would
-			 have to know to not free them)  */
-		for (i=0; i<ngaus; i++) {
-			(*pxgaus)[i]=xip [ngaus-1][i];
-			(*pxwgt )[i]=wgtp[ngaus-1][i];
-		}
-	}
-	else {
-
-		/*  calculate the Gauss points using recurrence relations  */
-		alpha=(double *) xmalloc(ngaus*sizeof(double));
-		beta =(double *) xmalloc(ngaus*sizeof(double));
-
-		/*  calculate the Legendre recurrence coefficients  */
-		alpha[0]=0.;
-		beta [0]=2.;
-
-		for (i=1; i<ngaus; i++) {
-			alpha[i]=0.;
-			beta [i]=1./(4.-1./(i*i));
-		}
-
-		/*  calculate the Gauss-Lobatto quadrature formula  */
-		for (i=0; i<ngaus-1; i++) {
-			pm1l=p0l;
-			p0l=p1l;
-			pm1r=p0r;
-			p0r=p1r;
-			p1l=(left -alpha[i])*p0l-beta[i]*pm1l;
-			p1r=(right-alpha[i])*p0r-beta[i]*pm1r;
-		}
-
-		/*  Normalize system to prevent underflow:
-			 [ p1l p0l ]{ a } = {left *p1l}
-			 [ p1r p0r ]{ b }   {right*p1r}
-			 dividing by p1l in the first equation and p1r in the second.  */
-
-		//		det=p1l*p0r-p1r*p0l;
-		det=p0r/p1r-p0l/p1l;
-		//		alpha[ngaus-1]=(left*p1l*p0r-right*p1r*p0l)/det;
-		//		beta [ngaus-1]=(right-left)*p1l*p1r/det;
-		alpha[ngaus-1]=(left *(p0r/p1r)-right*(p0l/p1l))/det;
-		beta [ngaus-1]=(right          -left           )/det;
-
-		/*  calculate the Gauss points  */
-		GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
-		xfree((void **)&beta );
-		xfree((void **)&alpha);
-	}
-
-}/*}}}1*/
-/*FUNCTION GaussRecur{{{1*/
-void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {
-	/*Gauss quadrature points from recursion coefficients.
-	 *
-	 *The routine uses the algorithm from the ORTHPOL routine GAUSS, which
-	 *finds the eigenvalues of a tridiagonal matrix.*/
-
-	/*Intermediaries*/
-	int i,j,k,l,m,ii,mml,iter;
-	double p,g,r,s,c,f,b;
-	double* work;
-
-	if (n==1){
-		zero[0]  =alpha[0];
-		weight[0]=beta[0];
-		return;
-	}
-
-	work=(double*)xmalloc(n*sizeof(double));
-
-	zero[0]  =alpha[0];
-	weight[0]=1.;
-	work[n-1]=0.;
-	for (i=1; i<=n-1; i++){
-		zero[i]=alpha[i];
-		work[i-1]=sqrt(beta[i]);
-		weight[i]=0;
-	}
-
-	for (l=0; l<=n-1; l++){
-		iter=0;
-		do {
-
-			/*  Look for a small subdiagonal element.  */
-			for (m=l; m<=n-1; m++) {
-				if (m == n-1) break;
-				if (fabs(work[m])
-							<= DBL_EPSILON*(fabs(zero[m])+fabs(zero[m+1])))
-				 break;
-			}
-			p=zero[l];
-			if (m==l) break;
-			++iter;
-
-			/*  Form shift.  */
-			g=(zero[l+1]-p)/(2.*work[l]);
-			r=sqrt(g*g+1.);
-			//			g=zero[m]-p+work[l]/(g+FortranSign(r,g));
-			g=zero[m]-p+work[l]/(g+(g>=0 ? fabs(r) : -fabs(r)));
-			s=1.;
-			c=1.;
-			p=0.;
-			mml=m-l;
-
-			/*  For i=m-1 step -1 until l do ...  */
-			for (ii=1; ii<=mml; ii++) {
-				i=m-ii;
-				f=s*work[i];
-				b=c*work[i];
-				if (fabs(f) >= fabs(g)) {
-					c=g/f;
-					r=sqrt(c*c+1.);
-					work[i+1]=f*r;
-					s=1./r;
-					c*=s;
-				}
-				else {
-					s=f/g;
-					r=sqrt(s*s+1.);
-					work[i+1]=g*r;
-					c=1./r;
-					s*=c;
-				}
-				g=zero[i+1]-p;
-				r=(zero[i]-g)*s+2.*c*b;
-				p=s*r;
-				zero[i+1]=g+p;
-				g=c*r-b;
-
-				/*  Form first component of vector.  */
-				f=weight[i+1];
-				weight[i+1]=s*weight[i]+c*f;
-				weight[i  ]=c*weight[i]-s*f;
-			}
-			zero[l]-=p;
-			work[l]=g;
-			work[m]=0.;
-		} while (iter < MAX_GAUS_ITER);
-		if (iter >= MAX_GAUS_ITER) {
-			xfree((void **)&work);
-			ISSMERROR("%s%i"," Max iterations exceeded for l=",MAX_GAUS_ITER);
-		}
-	}
-
-	/*  Order eigenvalues and eigenvectors.  */
-	for (i=0;i<n-1;i++) {
-		k=i;
-		p=zero[i];
-		for (j=i+1;j<n;j++){
-			if (zero[j] < p){
-				k=j;
-				p=zero[j];
-			}
-		}
-		if (k > i){
-			p=zero[i];
-			zero[i]=zero[k];
-			zero[k]=p;
-			p=weight[i];
-			weight[i]=weight[k];
-			weight[k]=p;
-		}
-	}
-	for (i=0;i<n;i++){
-		weight[i]=beta[0]*weight[i]*weight[i];
-	}
-
-	/*Cleanup*/
-	xfree((void **)&work);
-
-}/*}}}1*/
-
-/*Element Gauss points*/
-/*FUNCTION GaussQuad{{{1*/
-void GaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) { 
-	/*Gauss quadrature points for the quadrilaterial.*/
-
-	/*get the gauss points using the product of two line rules  */
-	GaussLegendre(pxgaus, pxwgt, nigaus);
-	GaussLegendre(pegaus, pewgt, njgaus);
-}/*}}}1*/
-/*FUNCTION GaussTria{{{1*/
-void GaussTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {
+/*FUNCTION GaussLegendreTria{{{1*/
+void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {
 	/*Gauss quadrature points for the triangle.
 
@@ -1379,5 +1139,5 @@
 		sizeof(wgt20)/sizeof(double)};
 
-	//	_printf_("GaussTria: iord=%d\n",iord);
+	//	_printf_("GaussLegendreTria: iord=%d\n",iord);
 
 	/*  check to see if Gauss points need to be calculated  */
@@ -1414,5 +1174,5 @@
 
 		/*  get the gauss points in each direction  */
-		GaussLegendre(&xgaus, &xwgt, nigaus);
+		GaussLegendreLinear(&xgaus, &xwgt, nigaus);
 
 		egaus=xgaus;
@@ -1439,5 +1199,5 @@
 	}
 
-	//	_printf_("GaussTria - ngaus=%d\n",*pngaus);
+	//	_printf_("GaussLegendreTria - ngaus=%d\n",*pngaus);
 	//	for (i=0; i<*pngaus; i++)
 	//		_printf_("i=%d: l1gaus=%f,l2gaus=%f,l3gaus=%f,wgt=%f\n",
@@ -1446,24 +1206,6 @@
 	return;
 }/*}}}1*/
-/*FUNCTION GaussHexa{{{1*/
-void GaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) {
-	/*Gauss quadrature points for the hexahedron.*/
-
-	/*  get the gauss points using the product of three line rules  */
-	GaussLegendre(pxgaus, pxwgt, nigaus);
-	GaussLegendre(pegaus, pewgt, njgaus);
-	GaussLegendre(pzgaus, pzwgt, nkgaus);
-}/*}}}1*/
-/*FUNCTION GaussPenta{{{1*/
-void GaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ) {
-	/*Gauss quadrature points for the pentahedron.*/
-
-	/*  get the gauss points using the product of triangular and line rules  */
-	GaussTria(pngaus, pl1, pl2, pl3, pwgt, iord);
-	GaussLegendre(pzgaus, pzwgt, nkgaus);
-	
-}/*}}}1*/
-/*FUNCTION GaussTetra{{{1*/
-void GaussTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {
+/*FUNCTION GaussLegendreTetra{{{1*/
+void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {
 	/* Gauss quadrature points for the tetrahedron.
 
@@ -1656,5 +1398,5 @@
 		sizeof(wgt6 )/sizeof(double)};
 
-	//	_printf_("GaussTetra: iord=%d\n",iord);
+	//	_printf_("GaussLegendreTetra: iord=%d\n",iord);
 
 	/*  check to see if Gauss points need to be calculated  */
@@ -1695,5 +1437,5 @@
 
 		/*  get the gauss points in each direction  */
-		GaussLegendre(&xgaus, &xwgt, nigaus);
+		GaussLegendreLinear(&xgaus, &xwgt, nigaus);
 
 		egaus=xgaus;
@@ -1729,94 +1471,266 @@
 	}
 }/*}}}1*/
-/*FUNCTION GaussSegment{{{1*/
-void GaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss){
-	/* GaussSegment: 
-	 * This routine computes gaussian points on a reference segment from -1 to 1.
-	 * The order p of integration depends on the number of gaussian points and can be computed from the polynomial degree p that needs to be integrated.
-	 * numgauss= (p+1) /2 */
-
-	/*WARNING: This is a copy of GaussLegendre: Merge the two!*/
-
-	double* segment_coord=NULL;
-	double* gauss_weights=NULL;
+/*FUNCTION GaussLobatto{{{1*/
+void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {
+	/*Gauss-Lobatto quadrature points.
+
+	  The recurrence coefficients for Legendre polynomials on (-1,1)
+	  are defined (from the ORTHPOL subroutine RECUR with ipoly=1) as:
+
+	  alpha(i)=0.
+	  beta (i)=1./(4.-1./(i-1)^2))
+
+	  and then modified for the Gauss-Lobatto quadrature rule on (-1,1)
+	  (from the ORTHPOL subroutine LOB).
+
+	  For degree p, the required number of Gauss-Lobatto points is
+	  n>=(p+1)/2+1 (one more than Gauss-Legendre).*/
+
+	int i;
+	double *alpha,*beta;
+	double left=-1.,right= 1.;
+	double p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;
+
+	/*p= 1, npoint= 1 (Gauss-Legendre)*/
+	static double wgt1[]={2.000000000000000};
+	static double xi1[]={0.000000000000000};
+
+	/*p= 1, npoint= 2*/
+	static double wgt2[]={1.000000000000000, 1.000000000000000};
+	static double xi2[]={-1.000000000000000, 1.000000000000000};
+
+	/*p= 3, npoint= 3*/
+	static double wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
+	static double xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
+
+	/*p= 5, npoint= 4*/
+	static double wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};
+	static double xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};
+
+	/*p= 7, npoint= 5*/
+	static double wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};
+	static double xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};
+
+	static double* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };
+	static double* xip [MAX_LINE_GLOB_PTS]={xi1  ,xi2  ,xi3  ,xi4  ,xi5  };
+
+	static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof(double),
+		sizeof(wgt2 )/sizeof(double),
+		sizeof(wgt3 )/sizeof(double),
+		sizeof(wgt4 )/sizeof(double),
+		sizeof(wgt5 )/sizeof(double)};
+
+	//	_printf_("Gauss-Lobatto recurrence coefficients ngaus=%d\n",ngaus);
+	*pxgaus = (double *) xmalloc(ngaus*sizeof(double));
+	*pxwgt  = (double *) xmalloc(ngaus*sizeof(double));
+
+	/*  check to see if Gauss points need to be calculated  */
+	if (ngaus <= MAX_LINE_GLOB_PTS) {
+
+		/*  copy the points from the static arrays (noting that the pointers
+			 could be returned directly, but then the calling function would
+			 have to know to not free them)  */
+		for (i=0; i<ngaus; i++) {
+			(*pxgaus)[i]=xip [ngaus-1][i];
+			(*pxwgt )[i]=wgtp[ngaus-1][i];
+		}
+	}
+	else {
+
+		/*  calculate the Gauss points using recurrence relations  */
+		alpha=(double *) xmalloc(ngaus*sizeof(double));
+		beta =(double *) xmalloc(ngaus*sizeof(double));
+
+		/*  calculate the Legendre recurrence coefficients  */
+		alpha[0]=0.;
+		beta [0]=2.;
+
+		for (i=1; i<ngaus; i++) {
+			alpha[i]=0.;
+			beta [i]=1./(4.-1./(i*i));
+		}
+
+		/*  calculate the Gauss-Lobatto quadrature formula  */
+		for (i=0; i<ngaus-1; i++) {
+			pm1l=p0l;
+			p0l=p1l;
+			pm1r=p0r;
+			p0r=p1r;
+			p1l=(left -alpha[i])*p0l-beta[i]*pm1l;
+			p1r=(right-alpha[i])*p0r-beta[i]*pm1r;
+		}
+
+		/*  Normalize system to prevent underflow:
+			 [ p1l p0l ]{ a } = {left *p1l}
+			 [ p1r p0r ]{ b }   {right*p1r}
+			 dividing by p1l in the first equation and p1r in the second.  */
+
+		//		det=p1l*p0r-p1r*p0l;
+		det=p0r/p1r-p0l/p1l;
+		//		alpha[ngaus-1]=(left*p1l*p0r-right*p1r*p0l)/det;
+		//		beta [ngaus-1]=(right-left)*p1l*p1r/det;
+		alpha[ngaus-1]=(left *(p0r/p1r)-right*(p0l/p1l))/det;
+		beta [ngaus-1]=(right          -left           )/det;
+
+		/*  calculate the Gauss points  */
+		GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
+		xfree((void **)&beta );
+		xfree((void **)&alpha);
+	}
+
+}/*}}}1*/
+/*FUNCTION GaussRecur{{{1*/
+void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {
+	/*Gauss quadrature points from recursion coefficients.
+	 *
+	 *The routine uses the algorithm from the ORTHPOL routine GAUSS, which
+	 *finds the eigenvalues of a tridiagonal matrix.*/
+
+	/*Intermediaries*/
+	int i,j,k,l,m,ii,mml,iter;
+	double p,g,r,s,c,f,b;
+	double* work;
+
+	if (n==1){
+		zero[0]  =alpha[0];
+		weight[0]=beta[0];
+		return;
+	}
+
+	work=(double*)xmalloc(n*sizeof(double));
+
+	zero[0]  =alpha[0];
+	weight[0]=1.;
+	work[n-1]=0.;
+	for (i=1; i<=n-1; i++){
+		zero[i]=alpha[i];
+		work[i-1]=sqrt(beta[i]);
+		weight[i]=0;
+	}
+
+	for (l=0; l<=n-1; l++){
+		iter=0;
+		do {
+
+			/*  Look for a small subdiagonal element.  */
+			for (m=l; m<=n-1; m++) {
+				if (m == n-1) break;
+				if (fabs(work[m])
+							<= DBL_EPSILON*(fabs(zero[m])+fabs(zero[m+1])))
+				 break;
+			}
+			p=zero[l];
+			if (m==l) break;
+			++iter;
+
+			/*  Form shift.  */
+			g=(zero[l+1]-p)/(2.*work[l]);
+			r=sqrt(g*g+1.);
+			//			g=zero[m]-p+work[l]/(g+FortranSign(r,g));
+			g=zero[m]-p+work[l]/(g+(g>=0 ? fabs(r) : -fabs(r)));
+			s=1.;
+			c=1.;
+			p=0.;
+			mml=m-l;
+
+			/*  For i=m-1 step -1 until l do ...  */
+			for (ii=1; ii<=mml; ii++) {
+				i=m-ii;
+				f=s*work[i];
+				b=c*work[i];
+				if (fabs(f) >= fabs(g)) {
+					c=g/f;
+					r=sqrt(c*c+1.);
+					work[i+1]=f*r;
+					s=1./r;
+					c*=s;
+				}
+				else {
+					s=f/g;
+					r=sqrt(s*s+1.);
+					work[i+1]=g*r;
+					c=1./r;
+					s*=c;
+				}
+				g=zero[i+1]-p;
+				r=(zero[i]-g)*s+2.*c*b;
+				p=s*r;
+				zero[i+1]=g+p;
+				g=c*r-b;
+
+				/*  Form first component of vector.  */
+				f=weight[i+1];
+				weight[i+1]=s*weight[i]+c*f;
+				weight[i  ]=c*weight[i]-s*f;
+			}
+			zero[l]-=p;
+			work[l]=g;
+			work[m]=0.;
+		} while (iter < MAX_GAUS_ITER);
+		if (iter >= MAX_GAUS_ITER) {
+			xfree((void **)&work);
+			ISSMERROR("%s%i"," Max iterations exceeded for l=",MAX_GAUS_ITER);
+		}
+	}
+
+	/*  Order eigenvalues and eigenvectors.  */
+	for (i=0;i<n-1;i++) {
+		k=i;
+		p=zero[i];
+		for (j=i+1;j<n;j++){
+			if (zero[j] < p){
+				k=j;
+				p=zero[j];
+			}
+		}
+		if (k > i){
+			p=zero[i];
+			zero[i]=zero[k];
+			zero[k]=p;
+			p=weight[i];
+			weight[i]=weight[k];
+			weight[k]=p;
+		}
+	}
+	for (i=0;i<n;i++){
+		weight[i]=beta[0]*weight[i]*weight[i];
+	}
+
+	/*Cleanup*/
+	xfree((void **)&work);
+
+}/*}}}1*/
+
+/*Element Gauss points TO BE REMOVED*/
+/*FUNCTION gaussQuad{{{1*/
+void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) { 
+	/*Gauss quadrature points for the quadrilaterial.*/
+
+	/*get the gauss points using the product of two line rules  */
+	GaussLegendreLinear(pxgaus, pxwgt, nigaus);
+	GaussLegendreLinear(pegaus, pewgt, njgaus);
+}/*}}}1*/
+/*FUNCTION gaussHexa{{{1*/
+void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) {
+	/*Gauss quadrature points for the hexahedron.*/
+
+	/*  get the gauss points using the product of three line rules  */
+	GaussLegendreLinear(pxgaus, pxwgt, nigaus);
+	GaussLegendreLinear(pegaus, pewgt, njgaus);
+	GaussLegendreLinear(pzgaus, pzwgt, nkgaus);
+}/*}}}1*/
+/*FUNCTION gaussPenta{{{1*/
+void gaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ) {
+	/*Gauss quadrature points for the pentahedron.*/
+
+	/*  get the gauss points using the product of triangular and line rules  */
+	GaussLegendreTria(pngaus, pl1, pl2, pl3, pwgt, iord);
+	GaussLegendreLinear(pzgaus, pzwgt, nkgaus);
 	
-	
-	if (num_gauss==1){
-		//order=1, num_gauss=1. Can integrate polynomials of degree 0 to 1
-
-		gauss_weights=(double*)xmalloc(sizeof(double));
-		gauss_weights[0]= 2.0;
-		
-		segment_coord=(double*)xmalloc(sizeof(double));
-		segment_coord[0]= 0.0;
-		
-	}
-	else if (num_gauss==2){
-		//order=2, num_gauss=3. Can integrate polynomials of degree 0 to 3
-		
-		gauss_weights=(double*)xmalloc(num_gauss*sizeof(double));
-		gauss_weights[0]=1.0;
-		gauss_weights[1]=1.0;
-	
-		segment_coord=(double*)xmalloc(num_gauss*sizeof(double));
-		segment_coord[0]=-0.577350269189626;
-		segment_coord[1]=0.577350269189626;
-	
-	}
-	else if (num_gauss==3){
-		//order=3, num_gauss=4. Can integrate polynomials of degree 0 to 5
-
-		gauss_weights=(double*)xmalloc(num_gauss*sizeof(double));
-		gauss_weights[0]=0.555555555555555556;
-		gauss_weights[1]=0.88888888888889;
-		gauss_weights[2]=0.555555555555555556;
-		
-		segment_coord=(double*)xmalloc(num_gauss*sizeof(double));
-		segment_coord[0]=-0.774596669241483;
-		segment_coord[1]=0.0;
-		segment_coord[2]=0.774596669241483;
-		
-	}
-	else if (num_gauss==4){
-		//order=3, num_gauss=4. Can integrate polynomials of degree 0 to 5
-
-		gauss_weights=(double*)xmalloc(num_gauss*sizeof(double));
-		gauss_weights[0]=0.347854845137454;
-		gauss_weights[1]=0.652145154862546;
-		gauss_weights[2]=0.652145154862546;
-		gauss_weights[3]=0.347854845137454;
-		
-		segment_coord=(double*)xmalloc(num_gauss*sizeof(double));
-		segment_coord[0]=-0.861136311594053;
-		segment_coord[1]=-0.339981043584856;
-		segment_coord[2]=0.339981043584856;
-		segment_coord[3]=0.861136311594053;
-		
-	}
-	else if (num_gauss==5){
-		//order=3, num_gauss=4. Can integrate polynomials of degree 0 to 5
-
-		gauss_weights=(double*)xmalloc(num_gauss*sizeof(double));
-		gauss_weights[0]=0.236926885056189;
-		gauss_weights[1]=0.478628670499366;
-		gauss_weights[2]=0.568888888888889;
-		gauss_weights[3]=0.478628670499366;
-		gauss_weights[4]=0.236926885056189;
-		
-		segment_coord=(double*)xmalloc(num_gauss*sizeof(double));
-		segment_coord[0]=-0.906179845938664;
-		segment_coord[1]=-0.538469310105683;
-		segment_coord[2]=0.0;
-		segment_coord[3]=0.538469310105683;
-		segment_coord[4]=0.906179845938664;
-		
-	}
-	else{
-		_printf_("GaussSegment error message: order not supported yet");
-	}
-
-	/*Assign output pointers: */
-	*pgauss_weights=gauss_weights;
-	*psegment_coord=segment_coord;
-	return;
 }/*}}}1*/
+/*FUNCTION gaussSegment{{{1*/
+void gaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss){
+	/* gaussSegment*/
+
+	GaussLegendreLinear(psegment_coord,pgauss_weights,num_gauss);
+}/*}}}1*/
Index: /issm/trunk/src/c/shared/Numerics/GaussPoints.h
===================================================================
--- /issm/trunk/src/c/shared/Numerics/GaussPoints.h	(revision 5623)
+++ /issm/trunk/src/c/shared/Numerics/GaussPoints.h	(revision 5624)
@@ -7,25 +7,18 @@
 
 #define    MAX_LINE_GAUS_PTS    4
-void GaussLegendre( double** pxgaus, double** pxwgt, int ngaus ); 
-
+void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus ); 
+#define    MAX_TRIA_SYM_ORD    20
+void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord );
+#define    MAX_TETRA_SYM_ORD    6
+void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord );
 #define    MAX_LINE_GLOB_PTS    5
 void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ); 
-
 #define MAX_GAUS_ITER   30
 void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta );
 
-void GaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );
-
-#define    MAX_TRIA_SYM_ORD    20
-void GaussTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord );
-
-void GaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );
-
-void GaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus );
-
-#define    MAX_TETRA_SYM_ORD    6
-void GaussTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord );
-
-void GaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss);
+void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );
+void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );
+void gaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus );
+void gaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss);
 
 #endif