Index: /issm/trunk-jpl/src/c/shared/Elements/Arrhenius.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/Arrhenius.cpp (revision 12474) +++ /issm/trunk-jpl/src/c/shared/Elements/Arrhenius.cpp (revision 12475) @@ -6,5 +6,5 @@ #include -double Arrhenius(double temperature,double depth,double n){ +IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n){ /*Use EISMINT Parameterization for the rheology: Payne2000 * @@ -25,9 +25,9 @@ /*Some physical constants (Payne2000)*/ - double beta=8.66*pow(10.,-4.); - double R=8.314; + IssmDouble beta=8.66*pow(10.,-4.); + IssmDouble R=8.314; /*Intermediaries*/ - double A,B,Tstar; + IssmDouble A,B,Tstar; /*convert temperature to absolute temperature*/ Index: /issm/trunk-jpl/src/c/shared/Elements/CoordinateSystemTransform.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/CoordinateSystemTransform.cpp (revision 12474) +++ /issm/trunk-jpl/src/c/shared/Elements/CoordinateSystemTransform.cpp (revision 12475) @@ -5,12 +5,12 @@ #include -void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int* cs_array){ +void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array){ int i,counter; int numdofs = 0; - double norm; - double *transform = NULL; - double *values = NULL; - double coord_system[3][3]; + IssmDouble norm; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; + IssmDouble coord_system[3][3]; /*Some checks in debugging mode*/ @@ -27,5 +27,5 @@ /*Allocate and initialize transform matrix*/ - transform=xNew(numdofs*numdofs); + transform=xNew(numdofs*numdofs); for(i=0;i -double Paterson(double temperature){ +#include "../../include/include.h" + +IssmDouble Paterson(IssmDouble temperature){ /*output: */ - double B; - double T; + IssmDouble B; + IssmDouble T; /*Switch to celsius from Kelvin: */ @@ -30,39 +32,39 @@ if(T<=-45.0){ - B=pow((double)10,(double)8)*(-0.000292866376675*pow(T+50,3)+ 0.011672640664130*pow(T+50,2) -0.325004442485481*(T+50)+ 6.524779401948101); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000292866376675*pow(T+50,3)+ 0.011672640664130*pow(T+50,2) -0.325004442485481*(T+50)+ 6.524779401948101); } else if((T>=-45.0) && (T<=-40.0)){ - B=pow((double)10,(double)8)*(-0.000292866376675*pow(T+45,3)+ 0.007279645014004*pow(T+45,2) -0.230243014094813*(T+45)+ 5.154964909039554); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000292866376675*pow(T+45,3)+ 0.007279645014004*pow(T+45,2) -0.230243014094813*(T+45)+ 5.154964909039554); } else if((T>=-40.0) && (T<=-35.0)){ - B=pow((double)10,(double)8)*(0.000072737147457*pow(T+40,3)+ 0.002886649363879*pow(T+40,2) -0.179411542205399*(T+40)+ 4.149132666831214); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(0.000072737147457*pow(T+40,3)+ 0.002886649363879*pow(T+40,2) -0.179411542205399*(T+40)+ 4.149132666831214); } else if((T>=-35.0) && (T<=-30.0)){ - B=pow((double)10,(double)8)*(-0.000086144770023*pow(T+35,3)+ 0.003977706575736*pow(T+35,2) -0.145089762507325*(T+35)+ 3.333333333333331); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000086144770023*pow(T+35,3)+ 0.003977706575736*pow(T+35,2) -0.145089762507325*(T+35)+ 3.333333333333331); } else if((T>=-30.0) && (T<=-25.0)){ - B=pow((double)10,(double)8)*(-0.000043984685769*pow(T+30,3)+ 0.002685535025386*pow(T+30,2) -0.111773554501713*(T+30)+ 2.696559088937191); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000043984685769*pow(T+30,3)+ 0.002685535025386*pow(T+30,2) -0.111773554501713*(T+30)+ 2.696559088937191); } else if((T>=-25.0) && (T<=-20.0)){ - B=pow((double)10,(double)8)*(-0.000029799523463*pow(T+25,3)+ 0.002025764738854*pow(T+25,2) -0.088217055680511*(T+25)+ 2.199331606342181); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000029799523463*pow(T+25,3)+ 0.002025764738854*pow(T+25,2) -0.088217055680511*(T+25)+ 2.199331606342181); } else if((T>=-20.0) && (T<=-15.0)){ - B=pow((double)10,(double)8)*(0.000136920904777*pow(T+20,3)+ 0.001578771886910*pow(T+20,2) -0.070194372551690*(T+20)+ 1.805165505978111); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(0.000136920904777*pow(T+20,3)+ 0.001578771886910*pow(T+20,2) -0.070194372551690*(T+20)+ 1.805165505978111); } else if((T>=-15.0) && (T<=-10.0)){ - B=pow((double)10,(double)8)*(-0.000899763781026*pow(T+15,3)+ 0.003632585458564*pow(T+15,2) -0.044137585824322*(T+15)+ 1.510778053489523); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.000899763781026*pow(T+15,3)+ 0.003632585458564*pow(T+15,2) -0.044137585824322*(T+15)+ 1.510778053489523); } else if((T>=-10.0) && (T<=-5.0)){ - B=pow((double)10,(double)8)*(0.001676964325070*pow(T+10,3)- 0.009863871256831*pow(T+10,2) -0.075294014815659*(T+10)+ 1.268434288203714); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(0.001676964325070*pow(T+10,3)- 0.009863871256831*pow(T+10,2) -0.075294014815659*(T+10)+ 1.268434288203714); } else if((T>=-5.0) && (T<=-2.0)){ - B=pow((double)10,(double)8)*(-0.003748937622487*pow(T+5,3)+0.015290593619213*pow(T+5,2) -0.048160403003748*(T+5)+ 0.854987973338348); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.003748937622487*pow(T+5,3)+0.015290593619213*pow(T+5,2) -0.048160403003748*(T+5)+ 0.854987973338348); } else if(T>=-2.0){ - B=pow((double)10,(double)8)*(-0.003748937622488*pow(T+2,3)-0.018449844983174*pow(T+2,2) -0.057638157095631*(T+2)+ 0.746900791092860); + B=pow((IssmPDouble)10,(IssmPDouble)8)*(-0.003748937622488*pow(T+2,3)-0.018449844983174*pow(T+2,2) -0.057638157095631*(T+2)+ 0.746900791092860); } /*B cannot be negative!*/ - if(B<0) B=pow((double)10,(double)6); + if(B<0) B=pow((IssmPDouble)10,(IssmPDouble)6); return B; Index: /issm/trunk-jpl/src/c/shared/Elements/TransformInvStiffnessMatrixCoord.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/TransformInvStiffnessMatrixCoord.cpp (revision 12474) +++ /issm/trunk-jpl/src/c/shared/Elements/TransformInvStiffnessMatrixCoord.cpp (revision 12475) @@ -23,6 +23,6 @@ int i,j; int numdofs = 0; - double *transform = NULL; - double *values = NULL; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; /*Get total number of dofs*/ @@ -36,5 +36,5 @@ /*Copy current stiffness matrix*/ - values=xNew(Ke->nrows*Ke->ncols); + values=xNew(Ke->nrows*Ke->ncols); for(i=0;inrows;i++) for(j=0;jncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j]; @@ -49,5 +49,5 @@ /*Free Matrix*/ - xDelete(transform); - xDelete(values); + xDelete(transform); + xDelete(values); } Index: /issm/trunk-jpl/src/c/shared/Elements/TransformLoadVectorCoord.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/TransformLoadVectorCoord.cpp (revision 12474) +++ /issm/trunk-jpl/src/c/shared/Elements/TransformLoadVectorCoord.cpp (revision 12475) @@ -22,6 +22,6 @@ int i,j; int numdofs = 0; - double *transform = NULL; - double *values = NULL; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; /*Get total number of dofs*/ @@ -35,5 +35,5 @@ /*Copy current load vector*/ - values=xNew(pe->nrows); + values=xNew(pe->nrows); for(i=0;inrows;i++) values[i]=pe->values[i]; @@ -47,5 +47,5 @@ /*Free Matrices*/ - xDelete(transform); - xDelete(values); + xDelete(transform); + xDelete(values); } Index: /issm/trunk-jpl/src/c/shared/Elements/TransformSolutionCoord.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/TransformSolutionCoord.cpp (revision 12474) +++ /issm/trunk-jpl/src/c/shared/Elements/TransformSolutionCoord.cpp (revision 12475) @@ -4,5 +4,5 @@ #include "./elements.h" -void TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int cs_enum){ +void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int cs_enum){ int* cs_array=NULL; @@ -19,10 +19,10 @@ } -void TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int* cs_array){ +void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int* cs_array){ int i,j; int numdofs = 0; - double *transform = NULL; - double *values = NULL; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; /*Get total number of dofs*/ @@ -36,5 +36,5 @@ /*Copy current solution vector*/ - values=xNew(numdofs); + values=xNew(numdofs); for(i=0;i(transform); - xDelete(values); + xDelete(transform); + xDelete(values); } Index: /issm/trunk-jpl/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp (revision 12474) +++ /issm/trunk-jpl/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp (revision 12475) @@ -23,6 +23,6 @@ int i,j; int numdofs = 0; - double *transform = NULL; - double *values = NULL; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; /*Get total number of dofs*/ @@ -36,5 +36,5 @@ /*Copy current stiffness matrix*/ - values=xNew(Ke->nrows*Ke->ncols); + values=xNew(Ke->nrows*Ke->ncols); for(i=0;inrows;i++) for(j=0;jncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j]; @@ -49,5 +49,5 @@ /*Free Matrix*/ - xDelete(transform); - xDelete(values); + xDelete(transform); + xDelete(values); } Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 12474) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 12475) @@ -11,12 +11,12 @@ class ElementVector; -double Paterson(double temperature); -double Arrhenius(double temperature,double depth,double n); -void GetVerticesCoordinates(double* xyz, Node** nodes, int numvertices); +IssmDouble Paterson(IssmDouble temperature); +IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n); +void GetVerticesCoordinates(IssmDouble* xyz, Node** nodes, int numvertices); int GetNumberOfDofs( Node** nodes,int numnodes,int setenum,int approximation_enum); int* GetLocalDofList( Node** nodes,int numnodes,int setenum,int approximation_enum); int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation_enum); #ifdef _HAVE_DIAGNOSTIC_ -void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int* cs_array); +void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array); void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum); void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array); @@ -25,9 +25,9 @@ void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int cs_enum); void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int* cs_array); -void TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int cs_enum); -void TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int* cs_array); +void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int cs_enum); +void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int* cs_array); #endif -inline void printarray(double* array,int lines,int cols=1){ +inline void printarray(IssmDouble* array,int lines,int cols=1){ printf("\n"); for(int i=0;i(idima*idimc); + dtemp=xNew(idima*idimc); MatrixMultiply(a,nrowa,ncola,itrna,b,nrowb,ncolb,itrnb,dtemp,0); MatrixMultiply(dtemp,idima,idimc,0,c,nrowc,ncolc,itrnc,d,iaddd); - xDelete(dtemp); + xDelete(dtemp); } @@ -79,9 +79,9 @@ else { - dtemp=xNew(idimb*idimd); + dtemp=xNew(idimb*idimd); MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,dtemp,0); MatrixMultiply(a,nrowa,ncola,itrna,dtemp,idimb,idimd,0,d,iaddd); - xDelete(dtemp); + xDelete(dtemp); } @@ -89,5 +89,5 @@ }/*}}}*/ /*FUNCTION MatrixMuliply {{{*/ -int MatrixMultiply( double* a, int nrowa, int ncola, int itrna, double* b, int nrowb, int ncolb, int itrnb, double* c, int iaddc ){ +int MatrixMultiply( IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int iaddc ){ /*MatrixMultiply Perform matrix multiplication a*b+c.*/ int noerr=1; @@ -158,5 +158,5 @@ }/*}}}*/ /*FUNCTION MatrixInverse {{{*/ -int MatrixInverse( double* a, int ndim, int nrow, double* b, int nvec, double* pdet ){ +int MatrixInverse( IssmDouble* a, int ndim, int nrow, IssmDouble* b, int nvec, IssmDouble* pdet ){ /* MatrixInverse Perform matrix inversion and linear equation solution. @@ -172,5 +172,5 @@ int i,j,k,ipt,jpt,irow,icol,ipiv,ncol; int *pivrc1,*pivrc2,*pindx; - double pivot,det,dtemp; + IssmDouble pivot,det,dtemp; if (!b && nvec) { @@ -333,6 +333,6 @@ return noerr; }/*}}}*/ -/*FUNCTION Matrix2x2Determinant(double* Adet,double* A) {{{*/ -void Matrix2x2Determinant(double* Adet,double* A){ +/*FUNCTION Matrix2x2Determinant(IssmDouble* Adet,IssmDouble* A) {{{*/ +void Matrix2x2Determinant(IssmDouble* Adet,IssmDouble* A){ /*Compute determinant of a 2x2 matrix*/ @@ -341,9 +341,9 @@ } /*}}}*/ -/*FUNCTION Matrix2x2Invert(double* Ainv,double* A) {{{*/ -void Matrix2x2Invert(double* Ainv,double* A){ +/*FUNCTION Matrix2x2Invert(IssmDouble* Ainv,IssmDouble* A) {{{*/ +void Matrix2x2Invert(IssmDouble* Ainv,IssmDouble* A){ /*Intermediaries*/ - double det; + IssmDouble det; /*Compute determinant*/ @@ -358,6 +358,6 @@ }/*}}}*/ -/*FUNCTION Matrix3x3Determinant(double* Adet,double* A) {{{*/ -void Matrix3x3Determinant(double* Adet,double* A){ +/*FUNCTION Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A) {{{*/ +void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A){ /*Compute determinant of a 3x3 matrix*/ @@ -366,9 +366,9 @@ } /*}}}*/ -/*FUNCTION Matrix3x3Invert(double* Ainv,double* A) {{{*/ -void Matrix3x3Invert(double* Ainv,double* A){ +/*FUNCTION Matrix3x3Invert(IssmDouble* Ainv,IssmDouble* A) {{{*/ +void Matrix3x3Invert(IssmDouble* Ainv,IssmDouble* A){ /*Intermediaries*/ - double det; + IssmDouble det; /*Compute determinant*/ @@ -388,6 +388,6 @@ }/*}}}*/ -/*FUNCTION MatrixTranspose(double* Adet,double* A) {{{*/ -void MatrixTranspose(double* tA,double* A, int nrows, int ncols){ +/*FUNCTION MatrixTranspose(IssmDouble* Adet,IssmDouble* A) {{{*/ +void MatrixTranspose(IssmDouble* tA,IssmDouble* A, int nrows, int ncols){ /*Transpose a n*m matrix*/ Index: /issm/trunk-jpl/src/c/shared/Matrix/matrix.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Matrix/matrix.h (revision 12474) +++ /issm/trunk-jpl/src/c/shared/Matrix/matrix.h (revision 12475) @@ -6,12 +6,14 @@ #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 ); -void Matrix2x2Invert(double* Ainv, double* A); -void Matrix2x2Determinant(double* Adet,double* A); -void Matrix3x3Invert(double* Ainv, double* A); -void Matrix3x3Determinant(double* Adet,double* A); -void MatrixTranspose(double* tA,double* A,int nrows, int ncols); +#include "../../include/include.h" + +int TripleMultiply( IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int nrowc, int ncolc, int itrnc, IssmDouble* d, int iaddd); +int MatrixMultiply( IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int iaddc ); +int MatrixInverse( IssmDouble* a, int ndim, int nrow, IssmDouble* b, int nvec, IssmDouble* pdet ); +void Matrix2x2Invert(IssmDouble* Ainv, IssmDouble* A); +void Matrix2x2Determinant(IssmDouble* Adet,IssmDouble* A); +void Matrix3x3Invert(IssmDouble* Ainv, IssmDouble* A); +void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A); +void MatrixTranspose(IssmDouble* tA,IssmDouble* A,int nrows, int ncols); #endif //ifndef _MATRIXUTILS_H_