Index: /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp	(revision 98)
+++ /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp	(revision 99)
@@ -41,9 +41,7 @@
 #undef __FUNCT__
 #define __FUNCT__ "GaussLegendre"
-int GaussLegendre( double** pxgaus,
-				   double** pxwgt,
-				   int ngaus )
-{
-	int i,ierr=0;
+void GaussLegendre( double** pxgaus, double** pxwgt, int ngaus ) {
+	
+	int i;
 	double *alpha,*beta;
 
@@ -124,20 +122,9 @@
 /*  calculate the Gauss points  */
 
-		if ((ierr=GaussRecur(*pxgaus,
-							 *pxwgt,
-							 ngaus,
-							 alpha,
-							 beta )))
-			_printf_("%s -- Error %d from GaussRecur for ngaus=%d.\n",
-					 __FUNCT__,ierr,ngaus);
-
+		GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
 		xfree((void **)&beta );
 		xfree((void **)&alpha);
 	}
 
-//	for (i=0; i<ngaus; i++)
-//		_printf_("i=%d: xgaus=%f,xwgt=%f\n",i,(*pxgaus)[i],(*pxwgt)[i]);
-
-	return(ierr);
 }
 
@@ -179,9 +166,7 @@
 #define __FUNCT__ "GaussLobatto"
 
-int GaussLobatto( double** pxgaus,
-				  double** pxwgt,
-				  int ngaus )
-{
-	int i,ierr=0;
+void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {
+
+	int i;
 	double *alpha,*beta;
 	double left=-1.,right= 1.;
@@ -297,20 +282,9 @@
 /*  calculate the Gauss points  */
 
-		if ((ierr=GaussRecur(*pxgaus,
-							 *pxwgt,
-							 ngaus,
-							 alpha,
-							 beta )))
-			_printf_("%s -- Error %d from GaussRecur for ngaus=%d.\n",
-					 __FUNCT__,ierr,ngaus);
-
+		GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
 		xfree((void **)&beta );
 		xfree((void **)&alpha);
 	}
 
-//	for (i=0; i<ngaus; i++)
-//		_printf_("i=%d: xgaus=%f,xwgt=%f\n",i,(*pxgaus)[i],(*pxwgt)[i]);
-
-	return(ierr);
 }
 
@@ -339,20 +313,13 @@
 #define __FUNCT__ "GaussRecur"
 
-int GaussRecur( double* zero, 
-				double* weight,
-				int n,
-				double* alpha, 
-				double* beta )
-{
+void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {
 	int i,j,k,l,m,ii,mml,iter;
 	double p,g,r,s,c,f,b;
 	double* work;
 
-//	_printf_("Gauss abscissas and weights n=%d\n",n);
-
 	if ( n == 1 ) {
 		zero[0]  =alpha[0];
 		weight[0]=beta[0];
-		return 1;
+		return ;
 	}
 
@@ -463,5 +430,4 @@
 	xfree((void **)&work);
 
-	return 1;
 }
 
@@ -492,39 +458,11 @@
 #define __FUNCT__ "GaussQuad"
 
-int GaussQuad( double** pxgaus,
-			   double** pxwgt,
-			   double** pegaus,
-			   double** pewgt,
-			   int nigaus,
-			   int njgaus )
-{ 
-	int i,ierr=0;
-
-//	_printf_("GaussQuad: nigaus=%d,njgaus=%d\n",
-//			 nigaus,njgaus);
+void GaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) { 
 
 /*  get the gauss points using the product of two line rules  */
 
-	if ((ierr=GaussLegendre(pxgaus,
-							pxwgt,
-							nigaus))) {
-		_printf_("%s -- Error %d from GaussLegendre for nigaus=%d.\n",
-				 __FUNCT__,ierr,nigaus);
-		return(ierr);
-	}
-//	for (i=0; i<nigaus; i++)
-//		_printf_("i=%d: xgaus=%f,xwgt=%f\n",i,(*pxgaus)[i],(*pxwgt)[i]);
-
-	if ((ierr=GaussLegendre(pegaus,
-							pewgt,
-							njgaus))) {
-		_printf_("%s -- Error %d from GaussLegendre for njgaus=%d.\n",
-				 __FUNCT__,ierr,njgaus);
-		return(ierr);
-	}
-//	for (i=0; i<njgaus; i++)
-//		_printf_("j=%d: egaus=%f,ewgt=%f\n",i,(*pegaus)[i],(*pewgt)[i]);
-
-	return 1;
+	GaussLegendre(pxgaus, pxwgt, nigaus);
+
+	GaussLegendre(pegaus, pewgt, njgaus);
 }
 
@@ -559,12 +497,7 @@
 #undef __FUNCT__
 #define __FUNCT__ "GaussTria"
-int GaussTria( int* pngaus,
-			   double** pl1,
-			   double** pl2,
-			   double** pl3,
-			   double** pwgt,
-			   int iord )
-{
-	int i,j,ipt,nigaus,ierr=0;
+void GaussTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {
+	
+	int i,j,ipt,nigaus;
 	double xi,eta;
 	double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt;
@@ -1659,11 +1592,5 @@
 /*  get the gauss points in each direction  */
 
-		if ((ierr=GaussLegendre(&xgaus,
-								&xwgt,
-								nigaus))) {
-			_printf_("%s -- Error %d from GaussLegendre for nigaus=%d.\n",
-					 __FUNCT__,ierr,nigaus);
-			return(ierr);
-		}
+		GaussLegendre(&xgaus, &xwgt, nigaus);
 
 		egaus=xgaus;
@@ -1697,5 +1624,5 @@
 //				 i,(*pl1 )[i],(*pl2 )[i],(*pl3 )[i],(*pwgt)[i]);
 
-	return 1;
+	return ;
 }
 
@@ -1730,52 +1657,13 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "GaussHexa"
-int GaussHexa( double** pxgaus,
-			   double** pxwgt,
-			   double** pegaus,
-			   double** pewgt,
-			   double** pzgaus,
-			   double** pzwgt,
-			   int nigaus,
-			   int njgaus,
-			   int nkgaus )
-{ 
-	int i,ierr=0;
-
-//	_printf_("GaussHexa: nigaus=%d,njgaus=%d,nkgaus=%d\n",
-//			 nigaus,njgaus,nkgaus);
-
-/*  get the gauss points using the product of three line rules  */
-
-	if ((ierr=GaussLegendre(pxgaus,
-							pxwgt,
-							nigaus))) {
-		_printf_("%s -- Error %d from GaussLegendre for nigaus=%d.\n",
-				 __FUNCT__,ierr,nigaus);
-		return(ierr);
-	}
-//	for (i=0; i<nigaus; i++)
-//		_printf_("i=%d: xgaus=%f,xwgt=%f\n",i,(*pxgaus)[i],(*pxwgt)[i]);
-
-	if ((ierr=GaussLegendre(pegaus,
-							pewgt,
-							njgaus))) {
-		_printf_("%s -- Error %d from GaussLegendre for njgaus=%d.\n",
-				 __FUNCT__,ierr,njgaus);
-		return(ierr);
-	}
-//	for (i=0; i<njgaus; i++)
-//		_printf_("j=%d: egaus=%f,ewgt=%f\n",i,(*pegaus)[i],(*pewgt)[i]);
-
-	if ((ierr=GaussLegendre(pzgaus,
-							pzwgt,
-							nkgaus))) {
-		_printf_("%s -- Error %d from GaussLegendre for nkgaus=%d.\n",
-				 __FUNCT__,ierr,nkgaus);
-		return(ierr);
-	}
-//	for (i=0; i<nkgaus; i++)
-//		_printf_("k=%d: zgaus=%f,zwgt=%f\n",i,(*pzgaus)[i],(*pzwgt)[i]);
-
-	return 1;
+void GaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) { 
+
+	/*  get the gauss points using the product of three line rules  */
+
+	GaussLegendre(pxgaus, pxwgt, nigaus);
+
+	GaussLegendre(pegaus, pewgt, njgaus);
+
+	GaussLegendre(pzgaus, pzwgt, nkgaus);
 }
 
@@ -1811,46 +1699,12 @@
 #define __FUNCT__ "GaussPenta"
 
-int GaussPenta( int* pngaus,
-				double** pl1,
-				double** pl2,
-				double** pl3,
-				double** pwgt,
-				double** pzgaus,
-				double** pzwgt,
-				int iord,
-				int nkgaus )
-{ 
-	int i,ierr=0;
-
-//	_printf_("GaussPenta: iord=%d,nkgaus=%d\n",
-//			 iord,nkgaus);
+void GaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ) {
 
 /*  get the gauss points using the product of triangular and line rules  */
 
-	if ((ierr=GaussTria(pngaus,
-						pl1,
-						pl2,
-						pl3,
-						pwgt,
-						iord))) {
-		_printf_("%s -- Error %d from GaussTria for iord=%d.\n",
-				 __FUNCT__,ierr,iord);
-		return(ierr);
-	}
-//	for (i=0; i<*pngaus; i++)
-//		_printf_("i=%d: l1gaus=%f,l2gaus=%f,l3gaus=%f,wgt=%f\n",
-//				 i,(*pl1 )[i],(*pl2 )[i],(*pl3 )[i],(*pwgt)[i]);
-
-	if ((ierr=GaussLegendre(pzgaus,
-							pzwgt,
-							nkgaus))) {
-		_printf_("%s -- Error %d from GaussLegendre for nkgaus=%d.\n",
-				 __FUNCT__,ierr,nkgaus);
-		return(ierr);
-	}
-//	for (i=0; i<nkgaus; i++)
-//		_printf_("k=%d: zgaus=%f,zwgt=%f\n",i,(*pzgaus)[i],(*pzwgt)[i]);
-
-	return 1;
+	GaussTria(pngaus, pl1, pl2, pl3, pwgt, iord);
+
+	GaussLegendre(pzgaus, pzwgt, nkgaus);
+	
 }
 
@@ -1891,13 +1745,6 @@
 #define __FUNCT__ "GaussTetra"
 
-int GaussTetra( int* pngaus,
-				double** pl1,
-				double** pl2,
-				double** pl3,
-				double** pl4,
-				double** pwgt,
-				int iord )
-{
-	int i,j,k,ipt,nigaus,ierr=0;
+void GaussTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {
+	int i,j,k,ipt,nigaus;
 	double xi,eta,zeta;
 	double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt;
@@ -2126,11 +1973,5 @@
 /*  get the gauss points in each direction  */
 
-		if ((ierr=GaussLegendre(&xgaus,
-								&xwgt,
-								nigaus))) {
-			_printf_("%s -- Error %d from GaussLegendre for nigaus=%d.\n",
-					 __FUNCT__,ierr,nigaus);
-			return(ierr);
-		}
+		GaussLegendre(&xgaus, &xwgt, nigaus);
 
 		egaus=xgaus;
@@ -2168,10 +2009,4 @@
 	}
 
-//	_printf_("GaussTetra - ngaus=%d\n",*pngaus);
-//	for (i=0; i<*pngaus; i++)
-//		_printf_("i=%d: l1gaus=%f,l2gaus=%f,l3gaus=%f,l4gaus=%f,wgt=%f\n",
-//				 i,(*pl1 )[i],(*pl2 )[i],(*pl3 )[i],(*pl4 )[i],(*pwgt)[i]);
-
-	return 1;
 }
 
Index: /issm/trunk/src/c/shared/Numerics/GaussPoints.h
===================================================================
--- /issm/trunk/src/c/shared/Numerics/GaussPoints.h	(revision 98)
+++ /issm/trunk/src/c/shared/Numerics/GaussPoints.h	(revision 99)
@@ -8,68 +8,24 @@
 
 #define    MAX_LINE_GAUS_PTS    4
-
-int GaussLegendre( double** pxgaus,
-				   double** pxwgt,
-				   int ngaus );
+void GaussLegendre( double** pxgaus, double** pxwgt, int ngaus ); 
 
 #define    MAX_LINE_GLOB_PTS    5
-
-int GaussLobatto( double** pxgaus,
-				  double** pxwgt,
-				  int ngaus );
+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 );
 
-int GaussRecur( double* zero,
-				double* weight,
-				int n,
-				double* alpha,
-				double* beta );
-
-int GaussQuad( double** pxgaus,
-			   double** pxwgt,
-			   double** pegaus,
-			   double** pewgt,
-			   int nigaus,
-			   int njgaus );
+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 );
 
-int 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 );
 
-int GaussHexa( double** pxgaus,
-			   double** pxwgt,
-			   double** pegaus,
-			   double** pewgt,
-			   double** pzgaus,
-			   double** pzwgt,
-			   int nigaus,
-			   int njgaus,
-			   int nkgaus );
-
-int GaussPenta( int* pngaus,
-				double** pl1,
-				double** pl2,
-				double** pl3,
-				double** pwgt,
-				double** pzgaus,
-				double** pzwgt,
-				int iord,
-				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 );
 
-int 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);
 
Index: /issm/trunk/src/c/shared/Numerics/cross.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/cross.cpp	(revision 99)
+++ /issm/trunk/src/c/shared/Numerics/cross.cpp	(revision 99)
@@ -0,0 +1,20 @@
+/*!\file:  cross.cpp
+ * \brief cross product for 2 vectors
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+void cross(double* result,double* vector1,double* vector2){
+
+	/*result,vector1 and vector2 are all assumed to be of size 3: */
+
+	result[0]=vector1[1]*vector2[2]-vector1[2]*vector2[1];
+	result[1]=vector1[2]*vector2[0]-vector1[0]*vector2[2];
+	result[2]=vector1[0]*vector2[1]-vector1[1]*vector2[0];
+
+}
+
Index: /issm/trunk/src/c/shared/Numerics/norm.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/norm.cpp	(revision 99)
+++ /issm/trunk/src/c/shared/Numerics/norm.cpp	(revision 99)
@@ -0,0 +1,17 @@
+/*!\file:  norm.cpp
+ * \brief vector norm
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "math.h"
+
+double norm(double* vector){
+
+	/*vector is assumed to be of size 3: */
+	return sqrt(pow(vector[0],2)+pow(vector[1],2)+pow(vector[2],2));
+}
Index: /issm/trunk/src/c/shared/Numerics/numerics.h
===================================================================
--- /issm/trunk/src/c/shared/Numerics/numerics.h	(revision 98)
+++ /issm/trunk/src/c/shared/Numerics/numerics.h	(revision 99)
@@ -13,4 +13,6 @@
 double OptFunc(double scalar, OptArgs* optargs);
 void BrentSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs);
+void cross(double* result,double* vector1,double* vector2);
+double norm(double* vector);
 
 #endif //ifndef _NUMERICS_H_