Changeset 4898


Ignore:
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)

Location:
issm/trunk/src/c
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r4885 r4898  
    40144014}
    40154015/*}}}*/
    4016 /*FUNCTION Penta::GetMatrixInvert {{{1*/
    4017 void Penta::GetMatrixInvert(double*  Ke_invert, double* Ke){
    4018         /*Inverse a 3 by 3 matrix only */
    4019 
    4020         double a,b,c,d,e,f,g,h,i;
    4021         double det;
    4022         int calculationdof=3;
    4023 
    4024         /*Take the matrix components: */
    4025         a=*(Ke+calculationdof*0+0);
    4026         b=*(Ke+calculationdof*0+1);
    4027         c=*(Ke+calculationdof*0+2);
    4028         d=*(Ke+calculationdof*1+0);
    4029         e=*(Ke+calculationdof*1+1);
    4030         f=*(Ke+calculationdof*1+2);
    4031         g=*(Ke+calculationdof*2+0);
    4032         h=*(Ke+calculationdof*2+1);
    4033         i=*(Ke+calculationdof*2+2);
    4034 
    4035         det=a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g);
    4036 
    4037         *(Ke_invert+calculationdof*0+0)=(e*i-f*h)/det;
    4038         *(Ke_invert+calculationdof*0+1)=(c*h-b*i)/det;
    4039         *(Ke_invert+calculationdof*0+2)=(b*f-c*e)/det;
    4040         *(Ke_invert+calculationdof*1+0)=(f*g-d*i)/det;
    4041         *(Ke_invert+calculationdof*1+1)=(a*i-c*g)/det;
    4042         *(Ke_invert+calculationdof*1+2)=(c*d-a*f)/det;
    4043         *(Ke_invert+calculationdof*2+0)=(d*h-e*g)/det;
    4044         *(Ke_invert+calculationdof*2+1)=(b*g-a*h)/det;
    4045         *(Ke_invert+calculationdof*2+2)=(a*e-b*d)/det;
    4046 
    4047 }
    4048 /*}}}*/
    40494016/*FUNCTION Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_coord) {{{1*/
    40504017void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_coord){
    … …  
    51865153
    51875154        /*Inverse the matrix corresponding to bubble part Kbb */
    5188         GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]);
     5155        Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);
    51895156
    51905157        /*Multiply matrices to create the reduce matrix Ke_reduced */
    … …  
    52305197
    52315198        /*Inverse the matrix corresponding to bubble part Kbb */
    5232         GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]);
     5199        Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);
    52335200
    52345201        /*Multiply matrices to create the reduce matrix Ke_reduced */

  • issm/trunk/src/c/objects/Elements/Penta.h

    r4885 r4898  
    140140                void      GetDofList1(int* doflist);
    141141                int     GetElementType(void);
    142                 void      GetMatrixInvert(double*  Ke_invert, double* Ke);
    143142                void      GetNodalFunctions(double* l1l6, double* gauss_coord);
    144143                void      GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord);

  • issm/trunk/src/c/objects/Elements/PentaRef.cpp

    r4885 r4898  
    698698        /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take
    699699         * the determinant of it: */
    700         const int NDOF3=3;
    701         double J[NDOF3][NDOF3];
    702 
     700        double J[3][3];
     701
     702        /*Get Jacobian*/
    703703        GetJacobian(&J[0][0],xyz_list,gauss);
    704704
    705         *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];
    706 
    707         if(*Jdet<0){
    708                 ISSMERROR("negative jacobian determinant!");
    709         }
     705        /*Get Determinant*/
     706        Matrix3x3Determinant(Jdet,&J[0][0]);
     707        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
     708
    710709}
    711710/*}}}*/
    712711/*FUNCTION PentaRef::GetJacobianInvert {{{1*/
    713 void PentaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
    714 
    715         double Jdet;
    716         const int NDOF3=3;
     712void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){
     713
     714        /*Jacobian*/
     715        double J[3][3];
    717716
    718717        /*Call Jacobian routine to get the jacobian:*/
    719         GetJacobian(Jinv, xyz_list, gauss);
     718        GetJacobian(&J[0][0], xyz_list, gauss);
    720719
    721720        /*Invert Jacobian matrix: */
    722         MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet);
     721        Matrix3x3Invert(Jinv,&J[0][0]);
    723722}
    724723/*}}}*/

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

    r4885 r4898  
    249249        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    250250         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    251 
    252         double x1,x2,x3,y1,y2,y3;
    253 
    254         x1=*(xyz_list+3*0+0);
    255         y1=*(xyz_list+3*0+1);
    256         x2=*(xyz_list+3*1+0);
    257         y2=*(xyz_list+3*1+1);
    258         x3=*(xyz_list+3*2+0);
    259         y3=*(xyz_list+3*2+1);
    260 
    261         *Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
    262 
    263         if(Jdet<0){
    264                 ISSMERROR("negative jacobian determinant!");
    265         }
     251        double J[2][2];
     252
     253        /*Get Jacobian*/
     254        GetJacobian(&J[0][0],xyz_list,gauss);
     255
     256        /*Get Determinant*/
     257        Matrix2x2Determinant(Jdet,&J[0][0]);
     258        if(Jdet<0) ISSMERROR("negative jacobian determinant!");
    266259
    267260}
    … …  
    296289void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
    297290
    298         double Jdet;
    299         const int NDOF2=2;
    300         const int numgrids=3;
     291        /*Jacobian*/
     292        double J[2][2];
    301293
    302294        /*Call Jacobian routine to get the jacobian:*/
    303         GetJacobian(Jinv, xyz_list, gauss);
     295        GetJacobian(&J[0][0], xyz_list, gauss);
    304296
    305297        /*Invert Jacobian matrix: */
    306         MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
     298        Matrix2x2Invert(Jinv,&J[0][0]);
    307299
    308300}

  • issm/trunk/src/c/shared/Matrix/MatrixUtils.cpp

    r4885 r4898  
    188188        }
    189189
     190        /*In debugging mode, check that we are not dealing with simple matrices*/
     191        ISSMASSERT(!(ndim==2 && nrow==2));
     192        ISSMASSERT(!(ndim==3 && nrow==3));
     193
    190194/*  initialize local variables and arrays  */
    191195
    … …  
    338342        return noerr;
    339343}/*}}}*/
     344/*FUNCTION Matrix2x2Determinant(double* Adet,double* A) {{{1*/
     345void Matrix2x2Determinant(double* Adet,double* A){
     346        /*Compute determinant of a 2x2 matrix*/
     347
     348        /*det = a*d - c*b*/
     349        *Adet= A[0]*A[3]-A[2]*A[1];
     350}
     351/*}}}*/
     352/*FUNCTION Matrix2x2Invert(double* Ainv,double* A) {{{1*/
     353void Matrix2x2Invert(double* Ainv,double* A){
     354
     355        /*Intermediaries*/
     356        double det;
     357
     358        /*Compute determinant*/
     359        Matrix2x2Determinant(&det,A);
     360        if (fabs(det) < DBL_EPSILON) ISSMERROR("Determinant smaller that machine epsilon");
     361
     362        /*Compute invert*/
     363        Ainv[0]=   A[3]/det; /* =  d/det */
     364        Ainv[1]= - A[1]/det; /* = -b/det */
     365        Ainv[2]= - A[2]/det; /* = -c/det */
     366        Ainv[3]=   A[0]/det; /* =  a/det */
     367
     368}/*}}}*/
     369/*FUNCTION Matrix3x3Determinant(double* Adet,double* A) {{{1*/
     370void Matrix3x3Determinant(double* Adet,double* A){
     371        /*Compute determinant of a 3x3 matrix*/
     372
     373        /*det = a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)*/
     374   *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];
     375}
     376/*}}}*/
     377/*FUNCTION Matrix3x3Invert(double* Ainv,double* A) {{{1*/
     378void Matrix3x3Invert(double* Ainv,double* A){
     379
     380        /*Intermediaries*/
     381        double det;
     382
     383        /*Compute determinant*/
     384        Matrix3x3Determinant(&det,A);
     385        if (fabs(det) < DBL_EPSILON) ISSMERROR("Determinant smaller that machine epsilon");
     386
     387        /*Compute invert*/
     388        Ainv[0]=(A[4]*A[8]-A[5]*A[7])/det; /* = (e*i-f*h)/det */
     389        Ainv[1]=(A[2]*A[7]-A[1]*A[8])/det; /* = (c*h-b*i)/det */
     390        Ainv[2]=(A[1]*A[5]-A[2]*A[4])/det; /* = (b*f-c*e)/det */
     391        Ainv[3]=(A[5]*A[6]-A[3]*A[8])/det; /* = (f*g-d*i)/det */
     392        Ainv[4]=(A[0]*A[8]-A[2]*A[6])/det; /* = (a*i-c*g)/det */
     393        Ainv[5]=(A[2]*A[3]-A[0]*A[5])/det; /* = (c*d-a*f)/det */
     394        Ainv[6]=(A[3]*A[7]-A[4]*A[6])/det; /* = (d*h-e*g)/det */
     395        Ainv[7]=(A[1]*A[6]-A[0]*A[7])/det; /* = (b*g-a*h)/det */
     396        Ainv[8]=(A[0]*A[4]-A[1]*A[3])/det; /* = (a*e-b*d)/det */
     397
     398}/*}}}*/

  • issm/trunk/src/c/shared/Matrix/matrix.h

    r3332 r4898  
    44
    55#ifndef _MATRIXUTILS_H_
    6 #define  _MATRIXUTILS_H_
     6#define _MATRIXUTILS_H_
    77
    8 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);
    9 int MatrixMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int iaddc );
    10 int MatrixInverse( double* a, int ndim, int nrow, double* b, int nvec, double* pdet );
    11        
     8int  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);
     9int  MatrixMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int iaddc );
     10int  MatrixInverse( double* a, int ndim, int nrow, double* b, int nvec, double* pdet );
     11void Matrix2x2Invert(double* Ainv, double* A);
     12void Matrix2x2Determinant(double* Adet,double* A);
     13void Matrix3x3Invert(double* Ainv, double* A);
     14void Matrix3x3Determinant(double* Adet,double* A);
    1215
    1316#endif //ifndef _MATRIXUTILS_H_
Note: See TracChangeset for help on using the changeset viewer.