Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4897)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4898)
@@ -4014,37 +4014,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetMatrixInvert {{{1*/
-void Penta::GetMatrixInvert(double*  Ke_invert, double* Ke){
-	/*Inverse a 3 by 3 matrix only */
-
-	double a,b,c,d,e,f,g,h,i;
-	double det;
-	int calculationdof=3;
-
-	/*Take the matrix components: */
-	a=*(Ke+calculationdof*0+0);
-	b=*(Ke+calculationdof*0+1);
-	c=*(Ke+calculationdof*0+2);
-	d=*(Ke+calculationdof*1+0);
-	e=*(Ke+calculationdof*1+1);
-	f=*(Ke+calculationdof*1+2);
-	g=*(Ke+calculationdof*2+0);
-	h=*(Ke+calculationdof*2+1);
-	i=*(Ke+calculationdof*2+2);
-
-	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;
-	*(Ke_invert+calculationdof*0+2)=(b*f-c*e)/det;
-	*(Ke_invert+calculationdof*1+0)=(f*g-d*i)/det;
-	*(Ke_invert+calculationdof*1+1)=(a*i-c*g)/det;
-	*(Ke_invert+calculationdof*1+2)=(c*d-a*f)/det;
-	*(Ke_invert+calculationdof*2+0)=(d*h-e*g)/det;
-	*(Ke_invert+calculationdof*2+1)=(b*g-a*h)/det;
-	*(Ke_invert+calculationdof*2+2)=(a*e-b*d)/det;
-
-}
-/*}}}*/
 /*FUNCTION Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_coord) {{{1*/
 void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_coord){
@@ -5186,5 +5153,5 @@
 
 	/*Inverse the matrix corresponding to bubble part Kbb */
-	GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]);
+	Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);
 
 	/*Multiply matrices to create the reduce matrix Ke_reduced */
@@ -5230,5 +5197,5 @@
 
 	/*Inverse the matrix corresponding to bubble part Kbb */
-	GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]);
+	Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);
 
 	/*Multiply matrices to create the reduce matrix Ke_reduced */
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4897)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4898)
@@ -140,5 +140,4 @@
 		void	  GetDofList1(int* doflist);
 		int     GetElementType(void);
-		void	  GetMatrixInvert(double*  Ke_invert, double* Ke);
 		void	  GetNodalFunctions(double* l1l6, double* gauss_coord);
 		void	  GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord);
Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 4897)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 4898)
@@ -698,27 +698,26 @@
 	/*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take 
 	 * the determinant of it: */
-	const int NDOF3=3;
-	double J[NDOF3][NDOF3];
-
+	double J[3][3];
+
+	/*Get Jacobian*/
 	GetJacobian(&J[0][0],xyz_list,gauss);
 
-	*Jdet= J[0][0]*J[1][1]*J[2][2]-J[0][0]*J[1][2]*J[2][1]-J[1][0]*J[0][1]*J[2][2]+J[1][0]*J[0][2]*J[2][1]+J[2][0]*J[0][1]*J[1][2]-J[2][0]*J[0][2]*J[1][1];
-
-	if(*Jdet<0){
-		ISSMERROR("negative jacobian determinant!");
-	}
+	/*Get Determinant*/
+	Matrix3x3Determinant(Jdet,&J[0][0]);
+	if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
+
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetJacobianInvert {{{1*/
-void PentaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
-
-	double Jdet;
-	const int NDOF3=3;
+void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){
+
+	/*Jacobian*/
+	double J[3][3];
 
 	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(Jinv, xyz_list, gauss);
+	GetJacobian(&J[0][0], xyz_list, gauss);
 
 	/*Invert Jacobian matrix: */
-	MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet);
+	Matrix3x3Invert(Jinv,&J[0][0]);
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 4897)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 4898)
@@ -249,19 +249,12 @@
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	double x1,x2,x3,y1,y2,y3;
-
-	x1=*(xyz_list+3*0+0);
-	y1=*(xyz_list+3*0+1);
-	x2=*(xyz_list+3*1+0);
-	y2=*(xyz_list+3*1+1);
-	x3=*(xyz_list+3*2+0);
-	y3=*(xyz_list+3*2+1);
-
-	*Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
-
-	if(Jdet<0){
-		ISSMERROR("negative jacobian determinant!");
-	}
+	double J[2][2];
+
+	/*Get Jacobian*/
+	GetJacobian(&J[0][0],xyz_list,gauss);
+
+	/*Get Determinant*/
+	Matrix2x2Determinant(Jdet,&J[0][0]);
+	if(Jdet<0) ISSMERROR("negative jacobian determinant!");
 
 }
@@ -296,13 +289,12 @@
 void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
 
-	double Jdet;
-	const int NDOF2=2;
-	const int numgrids=3;
+	/*Jacobian*/
+	double J[2][2];
 
 	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(Jinv, xyz_list, gauss);
+	GetJacobian(&J[0][0], xyz_list, gauss);
 
 	/*Invert Jacobian matrix: */
-	MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
+	Matrix2x2Invert(Jinv,&J[0][0]);
 
 }
Index: /issm/trunk/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- /issm/trunk/src/c/shared/Matrix/MatrixUtils.cpp	(revision 4897)
+++ /issm/trunk/src/c/shared/Matrix/MatrixUtils.cpp	(revision 4898)
@@ -188,4 +188,8 @@
 	}
 
+	/*In debugging mode, check that we are not dealing with simple matrices*/
+	ISSMASSERT(!(ndim==2 && nrow==2));
+	ISSMASSERT(!(ndim==3 && nrow==3));
+
 /*  initialize local variables and arrays  */
 
@@ -338,2 +342,57 @@
 	return noerr;
 }/*}}}*/
+/*FUNCTION Matrix2x2Determinant(double* Adet,double* A) {{{1*/
+void Matrix2x2Determinant(double* Adet,double* A){
+	/*Compute determinant of a 2x2 matrix*/
+
+	/*det = a*d - c*b*/
+	*Adet= A[0]*A[3]-A[2]*A[1];
+}
+/*}}}*/
+/*FUNCTION Matrix2x2Invert(double* Ainv,double* A) {{{1*/
+void Matrix2x2Invert(double* Ainv,double* A){
+
+	/*Intermediaries*/
+	double det;
+
+	/*Compute determinant*/
+	Matrix2x2Determinant(&det,A);
+	if (fabs(det) < DBL_EPSILON) ISSMERROR("Determinant smaller that machine epsilon");
+
+	/*Compute invert*/
+	Ainv[0]=   A[3]/det; /* =  d/det */
+	Ainv[1]= - A[1]/det; /* = -b/det */
+	Ainv[2]= - A[2]/det; /* = -c/det */
+	Ainv[3]=   A[0]/det; /* =  a/det */
+
+}/*}}}*/
+/*FUNCTION Matrix3x3Determinant(double* Adet,double* A) {{{1*/
+void Matrix3x3Determinant(double* Adet,double* A){
+	/*Compute determinant of a 3x3 matrix*/
+
+	/*det = a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)*/
+   *Adet= A[0]*A[4]*A[8]-A[0]*A[5]*A[7]-A[3]*A[1]*A[8]+A[3]*A[2]*A[7]+A[6]*A[1]*A[5]-A[6]*A[2]*A[4];
+}
+/*}}}*/
+/*FUNCTION Matrix3x3Invert(double* Ainv,double* A) {{{1*/
+void Matrix3x3Invert(double* Ainv,double* A){
+
+	/*Intermediaries*/
+	double det;
+
+	/*Compute determinant*/
+	Matrix3x3Determinant(&det,A);
+	if (fabs(det) < DBL_EPSILON) ISSMERROR("Determinant smaller that machine epsilon");
+
+	/*Compute invert*/
+	Ainv[0]=(A[4]*A[8]-A[5]*A[7])/det; /* = (e*i-f*h)/det */
+	Ainv[1]=(A[2]*A[7]-A[1]*A[8])/det; /* = (c*h-b*i)/det */
+	Ainv[2]=(A[1]*A[5]-A[2]*A[4])/det; /* = (b*f-c*e)/det */
+	Ainv[3]=(A[5]*A[6]-A[3]*A[8])/det; /* = (f*g-d*i)/det */
+	Ainv[4]=(A[0]*A[8]-A[2]*A[6])/det; /* = (a*i-c*g)/det */
+	Ainv[5]=(A[2]*A[3]-A[0]*A[5])/det; /* = (c*d-a*f)/det */
+	Ainv[6]=(A[3]*A[7]-A[4]*A[6])/det; /* = (d*h-e*g)/det */
+	Ainv[7]=(A[1]*A[6]-A[0]*A[7])/det; /* = (b*g-a*h)/det */
+	Ainv[8]=(A[0]*A[4]-A[1]*A[3])/det; /* = (a*e-b*d)/det */
+
+}/*}}}*/
Index: /issm/trunk/src/c/shared/Matrix/matrix.h
===================================================================
--- /issm/trunk/src/c/shared/Matrix/matrix.h	(revision 4897)
+++ /issm/trunk/src/c/shared/Matrix/matrix.h	(revision 4898)
@@ -4,10 +4,13 @@
 
 #ifndef _MATRIXUTILS_H_
-#define  _MATRIXUTILS_H_
+#define _MATRIXUTILS_H_
 
-int TripleMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int nrowc, int ncolc, int itrnc, double* d, int iaddd);
-int MatrixMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int iaddc );
-int MatrixInverse( double* a, int ndim, int nrow, double* b, int nvec, double* pdet );
-	
+int  TripleMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int nrowc, int ncolc, int itrnc, double* d, int iaddd);
+int  MatrixMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int iaddc );
+int  MatrixInverse( double* a, int ndim, int nrow, double* b, int nvec, double* pdet );
+void Matrix2x2Invert(double* Ainv, double* A);
+void Matrix2x2Determinant(double* Adet,double* A);
+void Matrix3x3Invert(double* Ainv, double* A);
+void Matrix3x3Determinant(double* Adet,double* A);
 
 #endif //ifndef _MATRIXUTILS_H_