Changeset 12424
- Timestamp:
- 06/15/12 14:36:38 (13 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/objects/Gauss/GaussTria.cpp
r12365 r12424 42 42 43 43 /*Intermediaties*/ 44 double *seg_coords = NULL;45 double *seg_weights = NULL;44 IssmPDouble *seg_coords = NULL; 45 IssmPDouble *seg_weights = NULL; 46 46 int i,index3; 47 47 … … 51 51 52 52 /*Allocate GaussTria fields*/ 53 coords1= (double*)xmalloc(numgauss*sizeof(double));54 coords2= (double*)xmalloc(numgauss*sizeof(double));55 coords3= (double*)xmalloc(numgauss*sizeof(double));56 weights= (double*)xmalloc(numgauss*sizeof(double));53 coords1=xNew<IssmDouble>(numgauss*sizeof(IssmDouble)); 54 coords2=xNew<IssmDouble>(numgauss*sizeof(IssmDouble)); 55 coords3=xNew<IssmDouble>(numgauss*sizeof(IssmDouble)); 56 weights=xNew<IssmDouble>(numgauss*sizeof(IssmDouble)); 57 57 58 58 /*Reverse index1 and 2 if necessary*/ … … 71 71 else if (index1==0 && index2==2){ 72 72 for(i=0;i<numgauss;i++) coords1[i]= 0.5*(1-seg_coords[i]); 73 for(i=0;i<numgauss;i++) coords2[i]= 0;73 for(i=0;i<numgauss;i++) coords2[i]= 0 ; 74 74 for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]); 75 75 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i]; … … 97 97 /*FUNCTION GaussTria::~GaussTria(){{{*/ 98 98 GaussTria::~GaussTria(){ 99 x free((void**)&weights);100 x free((void**)&coords1);101 x free((void**)&coords2);102 x free((void**)&coords3);99 xDelete<IssmDouble>(weights); 100 xDelete<IssmDouble>(coords1); 101 xDelete<IssmDouble>(coords2); 102 xDelete<IssmDouble>(coords3); 103 103 } 104 104 /*}}}*/ … … 199 199 /*}}}*/ 200 200 /*FUNCTION GaussTria::GaussFromCoords{{{*/ 201 void GaussTria::GaussFromCoords( double x,double y,double* xyz_list){201 void GaussTria::GaussFromCoords(IssmDouble x,IssmDouble y,IssmDouble* xyz_list){ 202 202 203 203 /*Intermediaries*/ 204 double area = 0;205 double x1,y1,x2,y2,x3,y3;204 IssmDouble area = 0; 205 IssmDouble x1,y1,x2,y2,x3,y3; 206 206 207 207 /*in debugging mode: check that the default constructor has been called*/ -
issm/trunk-jpl/src/c/objects/Gauss/GaussTria.h
r12365 r12424 15 15 private: 16 16 int numgauss; 17 double* weights;18 double* coords1;19 double* coords2;20 double* coords3;17 IssmDouble* weights; 18 IssmDouble* coords1; 19 IssmDouble* coords2; 20 IssmDouble* coords3; 21 21 22 22 public: 23 double weight;24 double coord1;25 double coord2;26 double coord3;23 IssmDouble weight; 24 IssmDouble coord1; 25 IssmDouble coord2; 26 IssmDouble coord3; 27 27 28 28 public: … … 38 38 int end(void); 39 39 void Echo(void); 40 void GaussFromCoords( double x1,double y1,double* xyz_list);40 void GaussFromCoords(IssmDouble x1,IssmDouble y1,IssmDouble* xyz_list); 41 41 void GaussPoint(int ig); 42 42 void GaussVertex(int iv); -
issm/trunk-jpl/src/c/objects/Loads/Friction.cpp
r12365 r12424 55 55 } 56 56 /*}}}*/ 57 /*FUNCTION Friction::GetAlpha2( double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/58 void Friction::GetAlpha2( double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){57 /*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/ 58 void Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){ 59 59 60 60 /*This routine calculates the basal friction coefficient … … 62 62 63 63 /*diverse: */ 64 double r,s;65 double drag_p, drag_q;66 double gravity,rho_ice,rho_water;67 double Neff;68 double thickness,bed;69 double vx,vy,vz,vmag;70 double drag_coefficient;71 double alpha2;64 IssmDouble r,s; 65 IssmDouble drag_p, drag_q; 66 IssmDouble gravity,rho_ice,rho_water; 67 IssmDouble Neff; 68 IssmDouble thickness,bed; 69 IssmDouble vx,vy,vz,vmag; 70 IssmDouble drag_coefficient; 71 IssmDouble alpha2; 72 72 73 73 /*Recover parameters: */ … … 119 119 } 120 120 /*}}}*/ 121 /*FUNCTION Friction::GetAlpha2( double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/122 void Friction::GetAlpha2( double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){121 /*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/ 122 void Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){ 123 123 124 124 /*This routine calculates the basal friction coefficient … … 126 126 127 127 /*diverse: */ 128 double r,s;129 double drag_p, drag_q;130 double gravity,rho_ice,rho_water;131 double Neff;132 double thickness,bed;133 double vx,vy,vz,vmag;134 double drag_coefficient;135 double alpha2;128 IssmDouble r,s; 129 IssmDouble drag_p, drag_q; 130 IssmDouble gravity,rho_ice,rho_water; 131 IssmDouble Neff; 132 IssmDouble thickness,bed; 133 IssmDouble vx,vy,vz,vmag; 134 IssmDouble drag_coefficient; 135 IssmDouble alpha2; 136 136 137 137 /*Recover parameters: */ … … 183 183 } 184 184 /*}}}*/ 185 /*FUNCTION Friction::GetAlphaComplement( double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/186 void Friction::GetAlphaComplement( double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){185 /*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/ 186 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){ 187 187 188 188 /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. … … 192 192 /*diverse: */ 193 193 int i; 194 double r,s;195 double vx,vy,vz,vmag;196 double drag_p,drag_q;197 double Neff;198 double drag_coefficient;199 double bed,thickness;200 double gravity,rho_ice,rho_water;201 double alpha_complement;194 IssmDouble r,s; 195 IssmDouble vx,vy,vz,vmag; 196 IssmDouble drag_p,drag_q; 197 IssmDouble Neff; 198 IssmDouble drag_coefficient; 199 IssmDouble bed,thickness; 200 IssmDouble gravity,rho_ice,rho_water; 201 IssmDouble alpha_complement; 202 202 203 203 /*Recover parameters: */ … … 249 249 } 250 250 /*}}}*/ 251 /*FUNCTION Friction::GetAlphaComplement( double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/252 void Friction::GetAlphaComplement( double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){251 /*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/ 252 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){ 253 253 254 254 /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. … … 258 258 /*diverse: */ 259 259 int i; 260 double r,s;261 double vx,vy,vz,vmag;262 double drag_p,drag_q;263 double Neff;264 double drag_coefficient;265 double bed,thickness;266 double gravity,rho_ice,rho_water;267 double alpha_complement;260 IssmDouble r,s; 261 IssmDouble vx,vy,vz,vmag; 262 IssmDouble drag_p,drag_q; 263 IssmDouble Neff; 264 IssmDouble drag_coefficient; 265 IssmDouble bed,thickness; 266 IssmDouble gravity,rho_ice,rho_water; 267 IssmDouble alpha_complement; 268 268 269 269 /*Recover parameters: */ … … 316 316 /*}}}*/ 317 317 /*FUNCTION Friction::GetInputValue{{{*/ 318 void Friction::GetInputValue( double* pvalue,GaussTria* gauss,int enum_type){318 void Friction::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type){ 319 319 320 320 Input* input=inputs->GetInput(enum_type); … … 325 325 /*}}}*/ 326 326 /*FUNCTION Friction::GetInputValue{{{*/ 327 void Friction::GetInputValue( double* pvalue,GaussPenta* gauss,int enum_type){327 void Friction::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type){ 328 328 329 329 Input* input=inputs->GetInput(enum_type); -
issm/trunk-jpl/src/c/objects/Loads/Friction.h
r12365 r12424 27 27 28 28 void Echo(void); 29 void GetAlpha2( double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);30 void GetAlpha2( double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);31 void GetAlphaComplement( double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum);32 void GetAlphaComplement( double* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);33 void GetInputValue( double* pvalue,GaussTria* gauss,int enum_type);34 void GetInputValue( double* pvalue,GaussPenta* gauss,int enum_type);29 void GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum); 30 void GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum); 31 void GetAlphaComplement(IssmDouble* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum); 32 void GetAlphaComplement(IssmDouble* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum); 33 void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type); 34 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type); 35 35 36 36 }; -
issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp
r12365 r12424 1 1 /* Gauss point structures and prototypes */ 2 2 3 #include "../../include/include.h" 3 4 #include "./GaussPoints.h" 4 5 #include "../Alloc/alloc.h" 5 #include "../../include/include.h"6 6 #include "../../io/io.h" 7 7 #include "../Exceptions/exceptions.h" … … 11 11 /*General Gauss points*/ 12 12 /*FUNCTION GaussLegendreLinear {{{*/ 13 void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus){13 void GaussLegendreLinear( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus){ 14 14 /* Gauss-Legendre quadrature points. 15 15 … … 25 25 /*Intermediaries*/ 26 26 int i; 27 double *alpha,*beta;27 IssmPDouble *alpha,*beta; 28 28 29 29 /*p= 1, npoint= 1*/ 30 static double wgt1[]={2.000000000000000};31 static double xi1[]={0.000000000000000};30 static IssmPDouble wgt1[]={2.000000000000000}; 31 static IssmPDouble xi1[]={0.000000000000000}; 32 32 33 33 /*p= 3, npoint= 2*/ 34 static double wgt2[]={1.000000000000000, 1.000000000000000};35 static double xi2[]={-0.577350269189626, 0.577350269189626};34 static IssmPDouble wgt2[]={1.000000000000000, 1.000000000000000}; 35 static IssmPDouble xi2[]={-0.577350269189626, 0.577350269189626}; 36 36 37 37 /*p= 5, npoint= 3*/ 38 static double wgt3[]={0.555555555555556, 0.888888888888889, 0.555555555555556};39 static double xi3[]={-0.774596669241483, 0.000000000000000, 0.774596669241483};38 static IssmPDouble wgt3[]={0.555555555555556, 0.888888888888889, 0.555555555555556}; 39 static IssmPDouble xi3[]={-0.774596669241483, 0.000000000000000, 0.774596669241483}; 40 40 41 41 /*p= 7, npoint= 4*/ 42 static double wgt4[]={0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454};43 static double xi4[]={-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053};44 45 static double* wgtp[MAX_LINE_GAUS_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 };46 static double* xip [MAX_LINE_GAUS_PTS]={xi1 ,xi2 ,xi3 ,xi4 };47 48 static int np[MAX_LINE_GAUS_PTS]={sizeof(wgt1 )/sizeof( double),49 sizeof(wgt2 )/sizeof( double),50 sizeof(wgt3 )/sizeof( double),51 sizeof(wgt4 )/sizeof( double)};42 static IssmPDouble wgt4[]={0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454}; 43 static IssmPDouble xi4[]={-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053}; 44 45 static IssmPDouble* wgtp[MAX_LINE_GAUS_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 }; 46 static IssmPDouble* xip [MAX_LINE_GAUS_PTS]={xi1 ,xi2 ,xi3 ,xi4 }; 47 48 static int np[MAX_LINE_GAUS_PTS]={sizeof(wgt1 )/sizeof(IssmPDouble), 49 sizeof(wgt2 )/sizeof(IssmPDouble), 50 sizeof(wgt3 )/sizeof(IssmPDouble), 51 sizeof(wgt4 )/sizeof(IssmPDouble)}; 52 52 53 53 // _printf_(true,"Gauss-Legendre recurrence coefficients ngaus=%d\n",ngaus); 54 *pxgaus = ( double *) xmalloc(ngaus*sizeof(double));55 *pxwgt = ( double *) xmalloc(ngaus*sizeof(double));54 *pxgaus = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble)); 55 *pxwgt = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble)); 56 56 57 57 /* check to see if Gauss points need to be calculated */ … … 70 70 71 71 /* calculate the Gauss points using recurrence relations */ 72 alpha=( double *) xmalloc(ngaus*sizeof(double));73 beta =( double *) xmalloc(ngaus*sizeof(double));72 alpha=(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble)); 73 beta =(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble)); 74 74 75 75 /* calculate the Legendre recurrence coefficients */ … … 89 89 }/*}}}*/ 90 90 /*FUNCTION GaussLegendreTria{{{*/ 91 void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {91 void GaussLegendreTria( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pwgt, int iord ) { 92 92 /*Gauss quadrature points for the triangle. 93 93 … … 98 98 /*Intermediaries*/ 99 99 int i,j,ipt,nigaus; 100 double xi,eta;101 double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt;100 IssmPDouble xi,eta; 101 IssmPDouble *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt; 102 102 103 103 /*Hardcoded Gauss points declaration*/ 104 104 /*p= 1, npoint= 1{{{*/ 105 static double wgt1[]={105 static IssmPDouble wgt1[]={ 106 106 1.732050807568877}; 107 static double l11[]={107 static IssmPDouble l11[]={ 108 108 0.333333333333333}; 109 static double l21[]={109 static IssmPDouble l21[]={ 110 110 0.333333333333333}; 111 static double l31[]={111 static IssmPDouble l31[]={ 112 112 0.333333333333333}; 113 113 /*}}}*/ 114 114 /*p= 2, npoint= 3 {{{*/ 115 static double wgt2[]={115 static IssmPDouble wgt2[]={ 116 116 0.577350269189625, 0.577350269189625, 0.577350269189625}; 117 static double l12[]={117 static IssmPDouble l12[]={ 118 118 0.666666666666667, 0.166666666666667, 0.166666666666667}; 119 static double l22[]={119 static IssmPDouble l22[]={ 120 120 0.166666666666667, 0.666666666666667, 0.166666666666667}; 121 static double l32[]={121 static IssmPDouble l32[]={ 122 122 0.166666666666667, 0.166666666666667, 0.666666666666667}; 123 123 /*}}}*/ 124 124 /*p= 3, npoint= 4 {{{*/ 125 static double wgt3[]={125 static IssmPDouble wgt3[]={ 126 126 -0.974278579257493, 0.902109795608790, 0.902109795608790, 127 127 0.902109795608790}; 128 static double l13[]={128 static IssmPDouble l13[]={ 129 129 0.333333333333333, 0.600000000000000, 0.200000000000000, 130 130 0.200000000000000}; 131 static double l23[]={131 static IssmPDouble l23[]={ 132 132 0.333333333333333, 0.200000000000000, 0.600000000000000, 133 133 0.200000000000000}; 134 static double l33[]={134 static IssmPDouble l33[]={ 135 135 0.333333333333333, 0.200000000000000, 0.200000000000000, 136 136 0.600000000000000}; 137 137 /*}}}*/ 138 138 /*p= 4, npoint= 6 {{{*/ 139 static double wgt4[]={139 static IssmPDouble wgt4[]={ 140 140 0.386908262797819, 0.386908262797819, 0.386908262797819, 141 141 0.190442006391807, 0.190442006391807, 0.190442006391807}; 142 static double l14[]={142 static IssmPDouble l14[]={ 143 143 0.108103018168070, 0.445948490915965, 0.445948490915965, 144 144 0.816847572980459, 0.091576213509771, 0.091576213509771}; 145 static double l24[]={145 static IssmPDouble l24[]={ 146 146 0.445948490915965, 0.108103018168070, 0.445948490915965, 147 147 0.091576213509771, 0.816847572980459, 0.091576213509771}; 148 static double l34[]={148 static IssmPDouble l34[]={ 149 149 0.445948490915965, 0.445948490915965, 0.108103018168070, 150 150 0.091576213509771, 0.091576213509771, 0.816847572980459}; 151 151 /*}}}*/ 152 152 /*p= 5, npoint= 7 {{{*/ 153 static double wgt5[]={153 static IssmPDouble wgt5[]={ 154 154 0.389711431702997, 0.229313399254729, 0.229313399254729, 155 155 0.229313399254729, 0.218133059367230, 0.218133059367230, 156 156 0.218133059367230}; 157 static double l15[]={157 static IssmPDouble l15[]={ 158 158 0.333333333333333, 0.059715871789770, 0.470142064105115, 159 159 0.470142064105115, 0.797426985353087, 0.101286507323456, 160 160 0.101286507323456}; 161 static double l25[]={161 static IssmPDouble l25[]={ 162 162 0.333333333333333, 0.470142064105115, 0.059715871789770, 163 163 0.470142064105115, 0.101286507323456, 0.797426985353087, 164 164 0.101286507323456}; 165 static double l35[]={165 static IssmPDouble l35[]={ 166 166 0.333333333333333, 0.470142064105115, 0.470142064105115, 167 167 0.059715871789770, 0.101286507323456, 0.101286507323456, … … 169 169 /*}}}*/ 170 170 /*p= 6, npoint=12 {{{*/ 171 static double wgt6[]={171 static IssmPDouble wgt6[]={ 172 172 0.202279763184836, 0.202279763184836, 0.202279763184836, 173 173 0.088065961139281, 0.088065961139281, 0.088065961139281, 174 174 0.143502272432755, 0.143502272432755, 0.143502272432755, 175 175 0.143502272432755, 0.143502272432755, 0.143502272432755}; 176 static double l16[]={176 static IssmPDouble l16[]={ 177 177 0.501426509658179, 0.249286745170910, 0.249286745170910, 178 178 0.873821971016996, 0.063089014491502, 0.063089014491502, 179 179 0.053145049844817, 0.053145049844817, 0.310352451033784, 180 180 0.636502499121399, 0.310352451033784, 0.636502499121399}; 181 static double l26[]={181 static IssmPDouble l26[]={ 182 182 0.249286745170910, 0.501426509658179, 0.249286745170910, 183 183 0.063089014491502, 0.873821971016996, 0.063089014491502, 184 184 0.310352451033784, 0.636502499121399, 0.053145049844817, 185 185 0.053145049844817, 0.636502499121399, 0.310352451033784}; 186 static double l36[]={186 static IssmPDouble l36[]={ 187 187 0.249286745170910, 0.249286745170910, 0.501426509658179, 188 188 0.063089014491502, 0.063089014491502, 0.873821971016996, … … 191 191 /*}}}*/ 192 192 /*p= 7, npoint=13 {{{*/ 193 static double wgt7[]={193 static IssmPDouble wgt7[]={ 194 194 -0.259062916308362, 0.304174548458604, 0.304174548458604, 195 195 0.304174548458604, 0.092400122517855, 0.092400122517855, … … 197 197 0.133564951824643, 0.133564951824643, 0.133564951824643, 198 198 0.133564951824643}; 199 static double l17[]={199 static IssmPDouble l17[]={ 200 200 0.333333333333333, 0.479308067841920, 0.260345966079040, 201 201 0.260345966079040, 0.869739794195568, 0.065130102902216, … … 203 203 0.312865496004874, 0.638444188569810, 0.312865496004874, 204 204 0.638444188569810}; 205 static double l27[]={205 static IssmPDouble l27[]={ 206 206 0.333333333333333, 0.260345966079040, 0.479308067841920, 207 207 0.260345966079040, 0.065130102902216, 0.869739794195568, … … 209 209 0.048690315425316, 0.048690315425316, 0.638444188569810, 210 210 0.312865496004874}; 211 static double l37[]={211 static IssmPDouble l37[]={ 212 212 0.333333333333333, 0.260345966079040, 0.260345966079040, 213 213 0.479308067841920, 0.065130102902216, 0.065130102902216, … … 217 217 /*}}}*/ 218 218 /*p= 8, npoint=16 {{{*/ 219 static double wgt8[]={219 static IssmPDouble wgt8[]={ 220 220 0.249961964823104, 0.164703541925695, 0.164703541925695, 221 221 0.164703541925695, 0.178777729989794, 0.178777729989794, … … 224 224 0.047164287656184, 0.047164287656184, 0.047164287656184, 225 225 0.047164287656184}; 226 static double l18[]={226 static IssmPDouble l18[]={ 227 227 0.333333333333333, 0.081414823414554, 0.459292588292723, 228 228 0.459292588292723, 0.658861384496480, 0.170569307751760, … … 231 231 0.263112829634638, 0.728492392955404, 0.263112829634638, 232 232 0.728492392955404}; 233 static double l28[]={233 static IssmPDouble l28[]={ 234 234 0.333333333333333, 0.459292588292723, 0.081414823414554, 235 235 0.459292588292723, 0.170569307751760, 0.658861384496480, … … 238 238 0.008394777409958, 0.008394777409958, 0.728492392955404, 239 239 0.263112829634638}; 240 static double l38[]={240 static IssmPDouble l38[]={ 241 241 0.333333333333333, 0.459292588292723, 0.459292588292723, 242 242 0.081414823414554, 0.170569307751760, 0.170569307751760, … … 247 247 /*}}}*/ 248 248 /*p= 9, npoint=19 {{{*/ 249 static double wgt9[]={249 static IssmPDouble wgt9[]={ 250 250 0.168244134395468, 0.054273292833345, 0.054273292833345, 251 251 0.054273292833345, 0.134801255248419, 0.134801255248419, … … 255 255 0.074969289332873, 0.074969289332873, 0.074969289332873, 256 256 0.074969289332873}; 257 static double l19[]={257 static IssmPDouble l19[]={ 258 258 0.333333333333333, 0.020634961602525, 0.489682519198738, 259 259 0.489682519198738, 0.125820817014127, 0.437089591492937, … … 263 263 0.221962989160766, 0.741198598784498, 0.221962989160766, 264 264 0.741198598784498}; 265 static double l29[]={265 static IssmPDouble l29[]={ 266 266 0.333333333333333, 0.489682519198738, 0.020634961602525, 267 267 0.489682519198738, 0.437089591492937, 0.125820817014127, … … 271 271 0.036838412054736, 0.036838412054736, 0.741198598784498, 272 272 0.221962989160766}; 273 static double l39[]={273 static IssmPDouble l39[]={ 274 274 0.333333333333333, 0.489682519198738, 0.489682519198738, 275 275 0.020634961602525, 0.437089591492937, 0.437089591492937, … … 281 281 /*}}}*/ 282 282 /*p=10, npoint=25 {{{*/ 283 static double wgt10[]={283 static IssmPDouble wgt10[]={ 284 284 0.157301373584232, 0.063611224790829, 0.063611224790829, 285 285 0.063611224790829, 0.078498377595183, 0.078498377595183, … … 291 291 0.016318805873179, 0.016318805873179, 0.016318805873179, 292 292 0.016318805873179}; 293 static double l110[]={293 static IssmPDouble l110[]={ 294 294 0.333333333333333, 0.028844733232685, 0.485577633383657, 295 295 0.485577633383657, 0.781036849029926, 0.109481575485037, … … 301 301 0.066803251012200, 0.923655933587500, 0.066803251012200, 302 302 0.923655933587500}; 303 static double l210[]={303 static IssmPDouble l210[]={ 304 304 0.333333333333333, 0.485577633383657, 0.028844733232685, 305 305 0.485577633383657, 0.109481575485037, 0.781036849029926, … … 311 311 0.009540815400299, 0.009540815400299, 0.923655933587500, 312 312 0.066803251012200}; 313 static double l310[]={313 static IssmPDouble l310[]={ 314 314 0.333333333333333, 0.485577633383657, 0.485577633383657, 315 315 0.028844733232685, 0.109481575485037, 0.109481575485037, … … 323 323 /*}}}*/ 324 324 /*p=11, npoint=27 {{{*/ 325 static double wgt11[]={325 static IssmPDouble wgt11[]={ 326 326 0.001605622060698, 0.001605622060698, 0.001605622060698, 327 327 0.133626914252765, 0.133626914252765, 0.133626914252765, … … 333 333 0.035866718600836, 0.035866718600836, 0.035866718600836, 334 334 0.035866718600836, 0.035866718600836, 0.035866718600836}; 335 static double l111[]={335 static IssmPDouble l111[]={ 336 336 -0.069222096541517, 0.534611048270758, 0.534611048270758, 337 337 0.202061394068290, 0.398969302965855, 0.398969302965855, … … 343 343 0.021022016536166, 0.021022016536166, 0.171488980304042, 344 344 0.807489003159792, 0.171488980304042, 0.807489003159792}; 345 static double l211[]={345 static IssmPDouble l211[]={ 346 346 0.534611048270758,-0.069222096541517, 0.534611048270758, 347 347 0.398969302965855, 0.202061394068290, 0.398969302965855, … … 353 353 0.171488980304042, 0.807489003159792, 0.021022016536166, 354 354 0.021022016536166, 0.807489003159792, 0.171488980304042}; 355 static double l311[]={355 static IssmPDouble l311[]={ 356 356 0.534611048270758, 0.534611048270758,-0.069222096541517, 357 357 0.398969302965855, 0.398969302965855, 0.202061394068290, … … 365 365 /*}}}*/ 366 366 /*p=12, npoint=33 {{{*/ 367 static double wgt12[]={367 static IssmPDouble wgt12[]={ 368 368 0.044567514407799, 0.044567514407799, 0.044567514407799, 369 369 0.075677707051848, 0.075677707051848, 0.075677707051848, … … 377 377 0.029992592075802, 0.029992592075802, 0.029992592075802, 378 378 0.029992592075802, 0.029992592075802, 0.029992592075802}; 379 static double l112[]={379 static IssmPDouble l112[]={ 380 380 0.023565220452390, 0.488217389773805, 0.488217389773805, 381 381 0.120551215411079, 0.439724392294460, 0.439724392294460, … … 389 389 0.025734050548330, 0.025734050548330, 0.116251915907597, 390 390 0.858014033544073, 0.116251915907597, 0.858014033544073}; 391 static double l212[]={391 static IssmPDouble l212[]={ 392 392 0.488217389773805, 0.023565220452390, 0.488217389773805, 393 393 0.439724392294460, 0.120551215411079, 0.439724392294460, … … 401 401 0.116251915907597, 0.858014033544073, 0.025734050548330, 402 402 0.025734050548330, 0.858014033544073, 0.116251915907597}; 403 static double l312[]={403 static IssmPDouble l312[]={ 404 404 0.488217389773805, 0.488217389773805, 0.023565220452390, 405 405 0.439724392294460, 0.439724392294460, 0.120551215411079, … … 415 415 /*}}}*/ 416 416 /* p=13, npoint=37 {{{*/ 417 static double wgt13[]={417 static IssmPDouble wgt13[]={ 418 418 0.090968907790622, 0.019537784619314, 0.019537784619314, 419 419 0.019537784619314, 0.054427130356344, 0.054427130356344, … … 429 429 0.026884523429480, 0.026884523429480, 0.026884523429480, 430 430 0.026884523429480}; 431 static double l113[]={431 static IssmPDouble l113[]={ 432 432 0.333333333333333, 0.009903630120591, 0.495048184939705, 433 433 0.495048184939705, 0.062566729780852, 0.468716635109574, … … 443 443 0.126357385491669, 0.851409537834241, 0.126357385491669, 444 444 0.851409537834241}; 445 static double l213[]={445 static IssmPDouble l213[]={ 446 446 0.333333333333333, 0.495048184939705, 0.009903630120591, 447 447 0.495048184939705, 0.468716635109574, 0.062566729780852, … … 457 457 0.022233076674090, 0.022233076674090, 0.851409537834241, 458 458 0.126357385491669}; 459 static double l313[]={459 static IssmPDouble l313[]={ 460 460 0.333333333333333, 0.495048184939705, 0.495048184939705, 461 461 0.009903630120591, 0.468716635109574, 0.468716635109574, … … 473 473 /*}}}*/ 474 474 /*p=14, npoint=42{{{*/ 475 static double wgt14[]={475 static IssmPDouble wgt14[]={ 476 476 0.037903474783419, 0.037903474783419, 0.037903474783419, 477 477 0.056791094234956, 0.056791094234956, 0.056791094234956, … … 488 488 0.008677970905831, 0.008677970905831, 0.008677970905831, 489 489 0.008677970905831, 0.008677970905831, 0.008677970905831}; 490 static double l114[]={490 static IssmPDouble l114[]={ 491 491 0.022072179275643, 0.488963910362179, 0.488963910362179, 492 492 0.164710561319092, 0.417644719340454, 0.417644719340454, … … 503 503 0.001268330932872, 0.001268330932872, 0.118974497696957, 504 504 0.879757171370171, 0.118974497696957, 0.879757171370171}; 505 static double l214[]={505 static IssmPDouble l214[]={ 506 506 0.488963910362179, 0.022072179275643, 0.488963910362179, 507 507 0.417644719340454, 0.164710561319092, 0.417644719340454, … … 518 518 0.118974497696957, 0.879757171370171, 0.001268330932872, 519 519 0.001268330932872, 0.879757171370171, 0.118974497696957}; 520 static double l314[]={520 static IssmPDouble l314[]={ 521 521 0.488963910362179, 0.488963910362179, 0.022072179275643, 522 522 0.417644719340454, 0.417644719340454, 0.164710561319092, … … 535 535 /*}}}*/ 536 536 /*p=15, npoint=48{{{*/ 537 static double wgt15[]={537 static IssmPDouble wgt15[]={ 538 538 0.003320126005206, 0.003320126005206, 0.003320126005206, 539 539 0.076641563419124, 0.076641563419124, 0.076641563419124, … … 552 552 0.013291658531346, 0.013291658531346, 0.013291658531346, 553 553 0.013291658531346, 0.013291658531346, 0.013291658531346}; 554 static double l115[]={554 static IssmPDouble l115[]={ 555 555 -0.013945833716486, 0.506972916858243, 0.506972916858243, 556 556 0.137187291433955, 0.431406354283023, 0.431406354283023, … … 569 569 0.012459809331199, 0.012459809331199, 0.103575616576386, 570 570 0.883964574092416, 0.103575616576386, 0.883964574092416}; 571 static double l215[]={571 static IssmPDouble l215[]={ 572 572 0.506972916858243,-0.013945833716486, 0.506972916858243, 573 573 0.431406354283023, 0.137187291433955, 0.431406354283023, … … 586 586 0.103575616576386, 0.883964574092416, 0.012459809331199, 587 587 0.012459809331199, 0.883964574092416, 0.103575616576386}; 588 static double l315[]={588 static IssmPDouble l315[]={ 589 589 0.506972916858243, 0.506972916858243,-0.013945833716486, 590 590 0.431406354283023, 0.431406354283023, 0.137187291433955, … … 605 605 /*}}}*/ 606 606 /*p=16, npoint=52 {{{*/ 607 static double wgt16[]={607 static IssmPDouble wgt16[]={ 608 608 0.081191089584902, 0.011095307165226, 0.011095307165226, 609 609 0.011095307165226, 0.072244353151393, 0.072244353151393, … … 624 624 0.011864642509229, 0.011864642509229, 0.011864642509229, 625 625 0.011864642509229}; 626 static double l116[]={626 static IssmPDouble l116[]={ 627 627 0.333333333333333, 0.005238916103123, 0.497380541948438, 628 628 0.497380541948438, 0.173061122901295, 0.413469438549352, … … 643 643 0.085283615682657, 0.900399064086661, 0.085283615682657, 644 644 0.900399064086661}; 645 static double l216[]={645 static IssmPDouble l216[]={ 646 646 0.333333333333333, 0.497380541948438, 0.005238916103123, 647 647 0.497380541948438, 0.413469438549352, 0.173061122901295, … … 662 662 0.014317320230681, 0.014317320230681, 0.900399064086661, 663 663 0.085283615682657}; 664 static double l316[]={664 static IssmPDouble l316[]={ 665 665 0.333333333333333, 0.497380541948438, 0.497380541948438, 666 666 0.005238916103123, 0.413469438549352, 0.413469438549352, … … 683 683 /*}}}*/ 684 684 /*p=17, npoint=61{{{*/ 685 static double wgt17[]={685 static IssmPDouble wgt17[]={ 686 686 0.057914928034477, 0.008822054327014, 0.008822054327014, 687 687 0.008822054327014, 0.025410682752829, 0.025410682752829, … … 705 705 0.011545213295771, 0.011545213295771, 0.011545213295771, 706 706 0.011545213295771}; 707 static double l117[]={707 static IssmPDouble l117[]={ 708 708 0.333333333333333, 0.005658918886452, 0.497170540556774, 709 709 0.497170540556774, 0.035647354750751, 0.482176322624625, … … 727 727 0.080711313679564, 0.904625504095608, 0.080711313679564, 728 728 0.904625504095608}; 729 static double l217[]={729 static IssmPDouble l217[]={ 730 730 0.333333333333333, 0.497170540556774, 0.005658918886452, 731 731 0.497170540556774, 0.482176322624625, 0.035647354750751, … … 749 749 0.014663182224828, 0.014663182224828, 0.904625504095608, 750 750 0.080711313679564}; 751 static double l317[]={751 static IssmPDouble l317[]={ 752 752 0.333333333333333, 0.497170540556774, 0.497170540556774, 753 753 0.005658918886452, 0.482176322624625, 0.482176322624625, … … 774 774 /* p=18, npoint=70 {{{*/ 775 775 776 static double wgt18[]={776 static IssmPDouble wgt18[]={ 777 777 0.053364381350150, 0.015713921277179, 0.015713921277179, 778 778 0.015713921277179, 0.032495554156279, 0.032495554156279, … … 799 799 0.000079999375178, 0.000079999375178, 0.000079999375178, 800 800 0.000079999375178}; 801 static double l118[]={801 static IssmPDouble l118[]={ 802 802 0.333333333333333, 0.013310382738157, 0.493344808630921, 803 803 0.493344808630921, 0.061578811516086, 0.469210594241957, … … 824 824 0.020874755282586, 1.014347260005363, 0.020874755282586, 825 825 1.014347260005363}; 826 static double l218[]={826 static IssmPDouble l218[]={ 827 827 0.333333333333333, 0.493344808630921, 0.013310382738157, 828 828 0.493344808630921, 0.469210594241957, 0.061578811516086, … … 849 849 -0.035222015287949,-0.035222015287949, 1.014347260005363, 850 850 0.020874755282586}; 851 static double l318[]={851 static IssmPDouble l318[]={ 852 852 0.333333333333333, 0.493344808630921, 0.493344808630921, 853 853 0.013310382738157, 0.469210594241957, 0.469210594241957, … … 877 877 /*p=19, npoint=73 {{{*/ 878 878 879 static double wgt19[]={879 static IssmPDouble wgt19[]={ 880 880 0.056995437856306, 0.017893352515055, 0.017893352515055, 881 881 0.017893352515055, 0.038775849701151, 0.038775849701151, … … 903 903 0.006581669842530, 0.006581669842530, 0.006581669842530, 904 904 0.006581669842530}; 905 static double l119[]={905 static IssmPDouble l119[]={ 906 906 0.333333333333333, 0.020780025853987, 0.489609987073006, 907 907 0.489609987073006, 0.090926214604215, 0.454536892697893, … … 929 929 0.065494628082938, 0.924344252620784, 0.065494628082938, 930 930 0.924344252620784}; 931 static double l219[]={931 static IssmPDouble l219[]={ 932 932 0.333333333333333, 0.489609987073006, 0.020780025853987, 933 933 0.489609987073006, 0.454536892697893, 0.090926214604215, … … 955 955 0.010161119296278, 0.010161119296278, 0.924344252620784, 956 956 0.065494628082938}; 957 static double l319[]={957 static IssmPDouble l319[]={ 958 958 0.333333333333333, 0.489609987073006, 0.489609987073006, 959 959 0.020780025853987, 0.454536892697893, 0.454536892697893, … … 983 983 /*}}}*/ 984 984 /*p=20, npoint=79 {{{*/ 985 static double wgt20[]={985 static IssmPDouble wgt20[]={ 986 986 0.057256499746719, 0.001501721280705, 0.001501721280705, 987 987 0.001501721280705, 0.020195803723819, 0.020195803723819, … … 1011 1011 0.006190192638113, 0.006190192638113, 0.006190192638113, 1012 1012 0.006190192638113}; 1013 static double l120[]={1013 static IssmPDouble l120[]={ 1014 1014 0.333333333333333,-0.001900928704400, 0.500950464352200, 1015 1015 0.500950464352200, 0.023574084130543, 0.488212957934729, … … 1039 1039 0.059696109149007, 0.929756171556853, 0.059696109149007, 1040 1040 0.929756171556853}; 1041 static double l220[]={1041 static IssmPDouble l220[]={ 1042 1042 0.333333333333333, 0.500950464352200,-0.001900928704400, 1043 1043 0.500950464352200, 0.488212957934729, 0.023574084130543, … … 1067 1067 0.010547719294141, 0.010547719294141, 0.929756171556853, 1068 1068 0.059696109149007}; 1069 static double l320[]={1069 static IssmPDouble l320[]={ 1070 1070 0.333333333333333, 0.500950464352200, 0.500950464352200, 1071 1071 -0.001900928704400, 0.488212957934729, 0.488212957934729, … … 1097 1097 /*}}}*/ 1098 1098 1099 static double* wgtp[MAX_TRIA_SYM_ORD]={1099 static IssmPDouble* wgtp[MAX_TRIA_SYM_ORD]={ 1100 1100 wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 , 1101 1101 wgt6 ,wgt7 ,wgt8 ,wgt9 ,wgt10, 1102 1102 wgt11,wgt12,wgt13,wgt14,wgt15, 1103 1103 wgt16,wgt17,wgt18,wgt19,wgt20}; 1104 static double* l1p [MAX_TRIA_SYM_ORD]={1104 static IssmPDouble* l1p [MAX_TRIA_SYM_ORD]={ 1105 1105 l11 ,l12 ,l13 ,l14 ,l15 , 1106 1106 l16 ,l17 ,l18 ,l19 ,l110 , 1107 1107 l111 ,l112 ,l113 ,l114 ,l115 , 1108 1108 l116 ,l117 ,l118 ,l119 ,l120 }; 1109 static double* l2p [MAX_TRIA_SYM_ORD]={1109 static IssmPDouble* l2p [MAX_TRIA_SYM_ORD]={ 1110 1110 l21 ,l22 ,l23 ,l24 ,l25 , 1111 1111 l26 ,l27 ,l28 ,l29 ,l210 , 1112 1112 l211 ,l212 ,l213 ,l214 ,l215 , 1113 1113 l216 ,l217 ,l218 ,l219 ,l220 }; 1114 static double* l3p [MAX_TRIA_SYM_ORD]={1114 static IssmPDouble* l3p [MAX_TRIA_SYM_ORD]={ 1115 1115 l31 ,l32 ,l33 ,l34 ,l35 , 1116 1116 l36 ,l37 ,l38 ,l39 ,l310 , … … 1118 1118 l316 ,l317 ,l318 ,l319 ,l320 }; 1119 1119 1120 static int np[MAX_TRIA_SYM_ORD]={sizeof(wgt1 )/sizeof( double),1121 sizeof(wgt2 )/sizeof( double),1122 sizeof(wgt3 )/sizeof( double),1123 sizeof(wgt4 )/sizeof( double),1124 sizeof(wgt5 )/sizeof( double),1125 sizeof(wgt6 )/sizeof( double),1126 sizeof(wgt7 )/sizeof( double),1127 sizeof(wgt8 )/sizeof( double),1128 sizeof(wgt9 )/sizeof( double),1129 sizeof(wgt10)/sizeof( double),1130 sizeof(wgt11)/sizeof( double),1131 sizeof(wgt12)/sizeof( double),1132 sizeof(wgt13)/sizeof( double),1133 sizeof(wgt14)/sizeof( double),1134 sizeof(wgt15)/sizeof( double),1135 sizeof(wgt16)/sizeof( double),1136 sizeof(wgt17)/sizeof( double),1137 sizeof(wgt18)/sizeof( double),1138 sizeof(wgt19)/sizeof( double),1139 sizeof(wgt20)/sizeof( double)};1120 static int np[MAX_TRIA_SYM_ORD]={sizeof(wgt1 )/sizeof(IssmPDouble), 1121 sizeof(wgt2 )/sizeof(IssmPDouble), 1122 sizeof(wgt3 )/sizeof(IssmPDouble), 1123 sizeof(wgt4 )/sizeof(IssmPDouble), 1124 sizeof(wgt5 )/sizeof(IssmPDouble), 1125 sizeof(wgt6 )/sizeof(IssmPDouble), 1126 sizeof(wgt7 )/sizeof(IssmPDouble), 1127 sizeof(wgt8 )/sizeof(IssmPDouble), 1128 sizeof(wgt9 )/sizeof(IssmPDouble), 1129 sizeof(wgt10)/sizeof(IssmPDouble), 1130 sizeof(wgt11)/sizeof(IssmPDouble), 1131 sizeof(wgt12)/sizeof(IssmPDouble), 1132 sizeof(wgt13)/sizeof(IssmPDouble), 1133 sizeof(wgt14)/sizeof(IssmPDouble), 1134 sizeof(wgt15)/sizeof(IssmPDouble), 1135 sizeof(wgt16)/sizeof(IssmPDouble), 1136 sizeof(wgt17)/sizeof(IssmPDouble), 1137 sizeof(wgt18)/sizeof(IssmPDouble), 1138 sizeof(wgt19)/sizeof(IssmPDouble), 1139 sizeof(wgt20)/sizeof(IssmPDouble)}; 1140 1140 1141 1141 // _printf_(true,"GaussLegendreTria: iord=%d\n",iord); … … 1150 1150 *pngaus=np[iord-1]; 1151 1151 1152 *pl1 = ( double *) xmalloc(*pngaus*sizeof(double));1153 *pl2 = ( double *) xmalloc(*pngaus*sizeof(double));1154 *pl3 = ( double *) xmalloc(*pngaus*sizeof(double));1155 *pwgt = ( double *) xmalloc(*pngaus*sizeof(double));1152 *pl1 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1153 *pl2 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1154 *pl3 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1155 *pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1156 1156 1157 1157 for (i=0; i<*pngaus; i++) { … … 1168 1168 *pngaus=nigaus*nigaus; 1169 1169 1170 *pl1 = ( double *) xmalloc(*pngaus*sizeof(double));1171 *pl2 = ( double *) xmalloc(*pngaus*sizeof(double));1172 *pl3 = ( double *) xmalloc(*pngaus*sizeof(double));1173 *pwgt = ( double *) xmalloc(*pngaus*sizeof(double));1170 *pl1 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1171 *pl2 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1172 *pl3 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1173 *pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1174 1174 1175 1175 /* get the gauss points in each direction */ … … 1207 1207 }/*}}}*/ 1208 1208 /*FUNCTION GaussLegendreTetra{{{*/ 1209 void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {1209 void GaussLegendreTetra( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pl4, IssmPDouble** pwgt, int iord ) { 1210 1210 /* Gauss quadrature points for the tetrahedron. 1211 1211 … … 1220 1220 /*Intermediaries*/ 1221 1221 int i,j,k,ipt,nigaus; 1222 double xi,eta,zeta;1223 double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt;1222 IssmPDouble xi,eta,zeta; 1223 IssmPDouble *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt; 1224 1224 1225 1225 /*Hardcoded Gauss points definition*/ 1226 1226 /*p= 1, npoint= 1 {{{*/ 1227 static double wgt1[]={1227 static IssmPDouble wgt1[]={ 1228 1228 1.000000000000000}; 1229 static double l11[]={1229 static IssmPDouble l11[]={ 1230 1230 0.250000000000000}; 1231 static double l21[]={1231 static IssmPDouble l21[]={ 1232 1232 0.250000000000000}; 1233 static double l31[]={1233 static IssmPDouble l31[]={ 1234 1234 0.250000000000000}; 1235 static double l41[]={1235 static IssmPDouble l41[]={ 1236 1236 0.250000000000000}; 1237 1237 /*}}}*/ 1238 1238 /*p= 2, npoint= 4 {{{*/ 1239 1239 1240 static double wgt2[]={1240 static IssmPDouble wgt2[]={ 1241 1241 0.250000000000000, 0.250000000000000, 0.250000000000000, 1242 1242 0.250000000000000}; 1243 static double l12[]={1243 static IssmPDouble l12[]={ 1244 1244 0.585410196624969, 0.138196601125011, 0.138196601125011, 1245 1245 0.138196601125011}; 1246 static double l22[]={1246 static IssmPDouble l22[]={ 1247 1247 0.138196601125011, 0.585410196624969, 0.138196601125011, 1248 1248 0.138196601125011}; 1249 static double l32[]={1249 static IssmPDouble l32[]={ 1250 1250 0.138196601125011, 0.138196601125011, 0.585410196624969, 1251 1251 0.138196601125011}; 1252 static double l42[]={1252 static IssmPDouble l42[]={ 1253 1253 0.138196601125011, 0.138196601125011, 0.138196601125011, 1254 1254 0.585410196624969}; 1255 1255 /*}}}*/ 1256 1256 /*p= 3, npoint= 5 {{{*/ 1257 static double wgt3[]={1257 static IssmPDouble wgt3[]={ 1258 1258 -0.800000000000000, 0.450000000000000, 0.450000000000000, 1259 1259 0.450000000000000, 0.450000000000000}; 1260 static double l13[]={1260 static IssmPDouble l13[]={ 1261 1261 0.250000000000000, 0.500000000000000, 0.166666666666667, 1262 1262 0.166666666666667, 0.166666666666667}; 1263 static double l23[]={1263 static IssmPDouble l23[]={ 1264 1264 0.250000000000000, 0.166666666666667, 0.500000000000000, 1265 1265 0.166666666666667, 0.166666666666667}; 1266 static double l33[]={1266 static IssmPDouble l33[]={ 1267 1267 0.250000000000000, 0.166666666666667, 0.166666666666667, 1268 1268 0.500000000000000, 0.166666666666667}; 1269 static double l43[]={1269 static IssmPDouble l43[]={ 1270 1270 0.250000000000000, 0.166666666666667, 0.166666666666667, 1271 1271 0.166666666666667, 0.500000000000000}; … … 1273 1273 /*p= 4, npoint=11 {{{*/ 1274 1274 1275 static double wgt4[]={1275 static IssmPDouble wgt4[]={ 1276 1276 -0.013155555555556, 0.007622222222222, 0.007622222222222, 1277 1277 0.007622222222222, 0.007622222222222, 0.024888888888889, 1278 1278 0.024888888888889, 0.024888888888889, 0.024888888888889, 1279 1279 0.024888888888889, 0.024888888888889}; 1280 static double l14[]={1280 static IssmPDouble l14[]={ 1281 1281 0.250000000000000, 0.785714285714286, 0.071428571428571, 1282 1282 0.071428571428571, 0.071428571428571, 0.399403576166799, 1283 1283 0.399403576166799, 0.399403576166799, 0.100596423833201, 1284 1284 0.100596423833201, 0.100596423833201}; 1285 static double l24[]={1285 static IssmPDouble l24[]={ 1286 1286 0.250000000000000, 0.071428571428571, 0.785714285714286, 1287 1287 0.071428571428571, 0.071428571428571, 0.399403576166799, 1288 1288 0.100596423833201, 0.100596423833201, 0.399403576166799, 1289 1289 0.399403576166799, 0.100596423833201}; 1290 static double l34[]={1290 static IssmPDouble l34[]={ 1291 1291 0.250000000000000, 0.071428571428571, 0.071428571428571, 1292 1292 0.785714285714286, 0.071428571428571, 0.100596423833201, 1293 1293 0.399403576166799, 0.100596423833201, 0.399403576166799, 1294 1294 0.100596423833201, 0.399403576166799}; 1295 static double l44[]={1295 static IssmPDouble l44[]={ 1296 1296 0.250000000000000, 0.071428571428571, 0.071428571428571, 1297 1297 0.071428571428571, 0.785714285714286, 0.100596423833201, … … 1301 1301 /*p= 5, npoint=15 {{{*/ 1302 1302 1303 static double wgt5[]={1303 static IssmPDouble wgt5[]={ 1304 1304 0.030283678097089, 0.006026785714286, 0.006026785714286, 1305 1305 0.006026785714286, 0.006026785714286, 0.011645249086029, … … 1307 1307 0.010949141561386, 0.010949141561386, 0.010949141561386, 1308 1308 0.010949141561386, 0.010949141561386, 0.010949141561386}; 1309 static double l15[]={1309 static IssmPDouble l15[]={ 1310 1310 0.250000000000000, 0.000000000000000, 0.333333333333333, 1311 1311 0.333333333333333, 0.333333333333333, 0.727272727272727, … … 1313 1313 0.066550153573664, 0.066550153573664, 0.066550153573664, 1314 1314 0.433449846426336, 0.433449846426336, 0.433449846426336}; 1315 static double l25[]={1315 static IssmPDouble l25[]={ 1316 1316 0.250000000000000, 0.333333333333333, 0.000000000000000, 1317 1317 0.333333333333333, 0.333333333333333, 0.090909090909091, … … 1319 1319 0.066550153573664, 0.433449846426336, 0.433449846426336, 1320 1320 0.066550153573664, 0.066550153573664, 0.433449846426336}; 1321 static double l35[]={1321 static IssmPDouble l35[]={ 1322 1322 0.250000000000000, 0.333333333333333, 0.333333333333333, 1323 1323 0.000000000000000, 0.333333333333333, 0.090909090909091, … … 1325 1325 0.433449846426336, 0.066550153573664, 0.433449846426336, 1326 1326 0.066550153573664, 0.433449846426336, 0.066550153573664}; 1327 static double l45[]={1327 static IssmPDouble l45[]={ 1328 1328 0.250000000000000, 0.333333333333333, 0.333333333333333, 1329 1329 0.333333333333333, 0.000000000000000, 0.090909090909091, … … 1334 1334 /*p= 6, npoint=24 {{{*/ 1335 1335 1336 static double wgt6[]={1336 static IssmPDouble wgt6[]={ 1337 1337 0.006653791709695, 0.006653791709695, 0.006653791709695, 1338 1338 0.006653791709695, 0.001679535175887, 0.001679535175887, … … 1343 1343 0.008035714285714, 0.008035714285714, 0.008035714285714, 1344 1344 0.008035714285714, 0.008035714285714, 0.008035714285714}; 1345 static double l16[]={1345 static IssmPDouble l16[]={ 1346 1346 0.356191386222545, 0.214602871259152, 0.214602871259152, 1347 1347 0.214602871259152, 0.877978124396166, 0.040673958534611, … … 1353 1353 0.269672331458316, 0.603005664791649, 0.269672331458316, 1354 1354 0.603005664791649, 0.269672331458316, 0.603005664791649}; 1355 static double l26[]={1355 static IssmPDouble l26[]={ 1356 1356 0.214602871259152, 0.356191386222545, 0.214602871259152, 1357 1357 0.214602871259152, 0.040673958534611, 0.877978124396166, … … 1363 1363 0.063661001875018, 0.063661001875018, 0.063661001875018, 1364 1364 0.063661001875018, 0.603005664791649, 0.269672331458316}; 1365 static double l36[]={1365 static IssmPDouble l36[]={ 1366 1366 0.214602871259152, 0.214602871259152, 0.356191386222545, 1367 1367 0.214602871259152, 0.040673958534611, 0.040673958534611, … … 1373 1373 0.063661001875018, 0.063661001875018, 0.603005664791649, 1374 1374 0.269672331458316, 0.063661001875018, 0.063661001875018}; 1375 static double l46[]={1375 static IssmPDouble l46[]={ 1376 1376 0.214602871259152, 0.214602871259152, 0.214602871259152, 1377 1377 0.356191386222545, 0.040673958534611, 0.040673958534611, … … 1385 1385 /*}}}*/ 1386 1386 1387 static double* wgtp[MAX_TETRA_SYM_ORD]={wgt1,wgt2,wgt3,wgt4,wgt5,wgt6};1388 static double* l1p [MAX_TETRA_SYM_ORD]={l11 ,l12 ,l13 ,l14 ,l15 ,l16 };1389 static double* l2p [MAX_TETRA_SYM_ORD]={l21 ,l22 ,l32 ,l24 ,l25 ,l26 };1390 static double* l3p [MAX_TETRA_SYM_ORD]={l31 ,l32 ,l33 ,l34 ,l35 ,l36 };1391 static double* l4p [MAX_TETRA_SYM_ORD]={l41 ,l42 ,l43 ,l44 ,l45 ,l46 };1392 1393 static int np[MAX_TETRA_SYM_ORD]={sizeof(wgt1 )/sizeof( double),1394 sizeof(wgt2 )/sizeof( double),1395 sizeof(wgt3 )/sizeof( double),1396 sizeof(wgt4 )/sizeof( double),1397 sizeof(wgt5 )/sizeof( double),1398 sizeof(wgt6 )/sizeof( double)};1387 static IssmPDouble* wgtp[MAX_TETRA_SYM_ORD]={wgt1,wgt2,wgt3,wgt4,wgt5,wgt6}; 1388 static IssmPDouble* l1p [MAX_TETRA_SYM_ORD]={l11 ,l12 ,l13 ,l14 ,l15 ,l16 }; 1389 static IssmPDouble* l2p [MAX_TETRA_SYM_ORD]={l21 ,l22 ,l32 ,l24 ,l25 ,l26 }; 1390 static IssmPDouble* l3p [MAX_TETRA_SYM_ORD]={l31 ,l32 ,l33 ,l34 ,l35 ,l36 }; 1391 static IssmPDouble* l4p [MAX_TETRA_SYM_ORD]={l41 ,l42 ,l43 ,l44 ,l45 ,l46 }; 1392 1393 static int np[MAX_TETRA_SYM_ORD]={sizeof(wgt1 )/sizeof(IssmPDouble), 1394 sizeof(wgt2 )/sizeof(IssmPDouble), 1395 sizeof(wgt3 )/sizeof(IssmPDouble), 1396 sizeof(wgt4 )/sizeof(IssmPDouble), 1397 sizeof(wgt5 )/sizeof(IssmPDouble), 1398 sizeof(wgt6 )/sizeof(IssmPDouble)}; 1399 1399 1400 1400 // _printf_(true,"GaussLegendreTetra: iord=%d\n",iord); … … 1410 1410 *pngaus=np[iord-1]; 1411 1411 1412 *pl1 = ( double *) xmalloc(*pngaus*sizeof(double));1413 *pl2 = ( double *) xmalloc(*pngaus*sizeof(double));1414 *pl3 = ( double *) xmalloc(*pngaus*sizeof(double));1415 *pl4 = ( double *) xmalloc(*pngaus*sizeof(double));1416 *pwgt = ( double *) xmalloc(*pngaus*sizeof(double));1412 *pl1 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1413 *pl2 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1414 *pl3 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1415 *pl4 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1416 *pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1417 1417 1418 1418 for (i=0; i<*pngaus; i++) { … … 1430 1430 *pngaus=nigaus*nigaus*nigaus; 1431 1431 1432 *pl1 = ( double *) xmalloc(*pngaus*sizeof(double));1433 *pl2 = ( double *) xmalloc(*pngaus*sizeof(double));1434 *pl3 = ( double *) xmalloc(*pngaus*sizeof(double));1435 *pl4 = ( double *) xmalloc(*pngaus*sizeof(double));1436 *pwgt = ( double *) xmalloc(*pngaus*sizeof(double));1432 *pl1 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1433 *pl2 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1434 *pl3 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1435 *pl4 = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1436 *pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble)); 1437 1437 1438 1438 /* get the gauss points in each direction */ … … 1472 1472 }/*}}}*/ 1473 1473 /*FUNCTION GaussLobatto{{{*/ 1474 void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {1474 void GaussLobatto( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus ) { 1475 1475 /*Gauss-Lobatto quadrature points. 1476 1476 … … 1488 1488 1489 1489 int i; 1490 double *alpha,*beta;1491 double left=-1.,right= 1.;1492 double p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;1490 IssmPDouble *alpha,*beta; 1491 IssmPDouble left=-1.,right= 1.; 1492 IssmPDouble p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det; 1493 1493 1494 1494 /*p= 1, npoint= 1 (Gauss-Legendre)*/ 1495 static double wgt1[]={2.000000000000000};1496 static double xi1[]={0.000000000000000};1495 static IssmPDouble wgt1[]={2.000000000000000}; 1496 static IssmPDouble xi1[]={0.000000000000000}; 1497 1497 1498 1498 /*p= 1, npoint= 2*/ 1499 static double wgt2[]={1.000000000000000, 1.000000000000000};1500 static double xi2[]={-1.000000000000000, 1.000000000000000};1499 static IssmPDouble wgt2[]={1.000000000000000, 1.000000000000000}; 1500 static IssmPDouble xi2[]={-1.000000000000000, 1.000000000000000}; 1501 1501 1502 1502 /*p= 3, npoint= 3*/ 1503 static double wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};1504 static double xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};1503 static IssmPDouble wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333}; 1504 static IssmPDouble xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000}; 1505 1505 1506 1506 /*p= 5, npoint= 4*/ 1507 static double wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};1508 static double xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};1507 static IssmPDouble wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667}; 1508 static IssmPDouble xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000}; 1509 1509 1510 1510 /*p= 7, npoint= 5*/ 1511 static double wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};1512 static double xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};1513 1514 static double* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };1515 static double* xip [MAX_LINE_GLOB_PTS]={xi1 ,xi2 ,xi3 ,xi4 ,xi5 };1516 1517 static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof( double),1518 sizeof(wgt2 )/sizeof( double),1519 sizeof(wgt3 )/sizeof( double),1520 sizeof(wgt4 )/sizeof( double),1521 sizeof(wgt5 )/sizeof( double)};1511 static IssmPDouble wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000}; 1512 static IssmPDouble xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000}; 1513 1514 static IssmPDouble* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 }; 1515 static IssmPDouble* xip [MAX_LINE_GLOB_PTS]={xi1 ,xi2 ,xi3 ,xi4 ,xi5 }; 1516 1517 static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof(IssmPDouble), 1518 sizeof(wgt2 )/sizeof(IssmPDouble), 1519 sizeof(wgt3 )/sizeof(IssmPDouble), 1520 sizeof(wgt4 )/sizeof(IssmPDouble), 1521 sizeof(wgt5 )/sizeof(IssmPDouble)}; 1522 1522 1523 1523 // _printf_(true,"Gauss-Lobatto recurrence coefficients ngaus=%d\n",ngaus); 1524 *pxgaus = ( double *) xmalloc(ngaus*sizeof(double));1525 *pxwgt = ( double *) xmalloc(ngaus*sizeof(double));1524 *pxgaus = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble)); 1525 *pxwgt = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble)); 1526 1526 1527 1527 /* check to see if Gauss points need to be calculated */ … … 1539 1539 1540 1540 /* calculate the Gauss points using recurrence relations */ 1541 alpha=( double *) xmalloc(ngaus*sizeof(double));1542 beta =( double *) xmalloc(ngaus*sizeof(double));1541 alpha=(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble)); 1542 beta =(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble)); 1543 1543 1544 1544 /* calculate the Legendre recurrence coefficients */ … … 1581 1581 }/*}}}*/ 1582 1582 /*FUNCTION GaussRecur{{{*/ 1583 void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {1583 void GaussRecur( IssmPDouble* zero, IssmPDouble* weight, int n, IssmPDouble* alpha, IssmPDouble* beta ) { 1584 1584 /*Gauss quadrature points from recursion coefficients. 1585 1585 * … … 1589 1589 /*Intermediaries*/ 1590 1590 int i,j,k,l,m,ii,mml,iter; 1591 double p,g,r,s,c,f,b;1592 double* work;1591 IssmPDouble p,g,r,s,c,f,b; 1592 IssmPDouble* work; 1593 1593 1594 1594 if (n==1){ … … 1598 1598 } 1599 1599 1600 work=( double*)xmalloc(n*sizeof(double));1600 work=(IssmPDouble*)xmalloc(n*sizeof(IssmPDouble)); 1601 1601 1602 1602 zero[0] =alpha[0]; … … 1704 1704 /*Element Gauss points TO BE REMOVED*/ 1705 1705 /*FUNCTION gaussQuad{{{*/ 1706 void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) {1706 void gaussQuad( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus ) { 1707 1707 /*Gauss quadrature points for the quadrilaterial.*/ 1708 1708 … … 1712 1712 }/*}}}*/ 1713 1713 /*FUNCTION gaussHexa{{{*/ 1714 void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) {1714 void gaussHexa( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus ) { 1715 1715 /*Gauss quadrature points for the hexahedron.*/ 1716 1716 -
issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.h
r5808 r12424 7 7 8 8 #define MAX_LINE_GAUS_PTS 4 9 void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus );9 void GaussLegendreLinear( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus ); 10 10 #define MAX_TRIA_SYM_ORD 20 11 void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord );11 void GaussLegendreTria( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pwgt, int iord ); 12 12 #define MAX_TETRA_SYM_ORD 6 13 void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord );13 void GaussLegendreTetra( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pl4, IssmPDouble** pwgt, int iord ); 14 14 #define MAX_LINE_GLOB_PTS 5 15 void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus );15 void GaussLobatto( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus ); 16 16 #define MAX_GAUS_ITER 30 17 void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta );17 void GaussRecur( IssmPDouble* zero, IssmPDouble* weight, int n, IssmPDouble* alpha, IssmPDouble* beta ); 18 18 19 void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );20 void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );19 void gaussQuad( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus ); 20 void gaussHexa( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus ); 21 21 22 22 #endif
Note:
See TracChangeset
for help on using the changeset viewer.