Changeset 4898 for issm/trunk/src/c/objects/Elements/TriaRef.cpp

Timestamp:

07/30/10 09:00:31 (15 years ago)

Author:

Mathieu Morlighem

Message:

Use Matrix utils to compute determinant and matrix invert. Simple routines for 2x2 and 3x3 matrices have been added (speeds up)

File:

• issm/trunk/src/c/objects/Elements/TriaRef.cpp (modified) (2 diffs)

issm/trunk/src/c/objects/Elements/TriaRef.cpp

@@ -249 +249 @@
         /*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 +289 @@
 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]);
 
 }