Ice Sheet System Model  4.18
Code documentation
Functions
GaussPoints.cpp File Reference
#include <math.h>
#include <float.h>
#include "./GaussPoints.h"
#include "../MemOps/MemOps.h"
#include "../Exceptions/exceptions.h"
#include "../Numerics/constants.h"

Go to the source code of this file.

Functions

void GaussLegendreLinear (IssmPDouble **pxgaus, IssmPDouble **pxwgt, int ngaus)
 
void GaussLegendreTria (int *pngaus, IssmDouble **pl1, IssmDouble **pl2, IssmDouble **pl3, IssmDouble **pwgt, int iord)
 
void GaussLegendreTetra (int *pngaus, IssmDouble **pl1, IssmDouble **pl2, IssmDouble **pl3, IssmDouble **pl4, IssmDouble **pwgt, int iord)
 
void GaussLobatto (IssmPDouble **pxgaus, IssmPDouble **pxwgt, int ngaus)
 
void GaussRecur (IssmPDouble *zero, IssmPDouble *weight, int n, IssmPDouble *alpha, IssmPDouble *beta)
 

Function Documentation

◆ GaussLegendreLinear()

void GaussLegendreLinear ( IssmPDouble **  pxgaus,
IssmPDouble **  pxwgt,
int  ngaus 
)

Definition at line 11 of file GaussPoints.cpp.

11  {/*{{{*/
12  /* Gauss-Legendre quadrature points.
13 
14  The recurrence coefficients for Legendre polynomials on (-1,1)
15  are defined (from the ORTHPOL subroutine RECUR with ipoly=1) as:
16 
17  alpha(i)=0.
18  beta (i)=1./(4.-1./(i-1)^2))
19 
20  For degree p, the required number of Gauss-Legendre points is
21  n>=(p+1)/2.*/
22  _assert_(ngaus>0);
23 
24  /*Intermediaries*/
25  int i;
26  IssmPDouble *alpha,*beta;
27 
28  /*p= 1, npoint= 1*/
29  static IssmPDouble wgt1[]={2.000000000000000};
30  static IssmPDouble xi1[]={0.000000000000000};
31 
32  /*p= 3, npoint= 2*/
33  static IssmPDouble wgt2[]={1.000000000000000, 1.000000000000000};
34  static IssmPDouble xi2[]={-0.577350269189626, 0.577350269189626};
35 
36  /*p= 5, npoint= 3*/
37  static IssmPDouble wgt3[]={0.555555555555556, 0.888888888888889, 0.555555555555556};
38  static IssmPDouble xi3[]={-0.774596669241483, 0.000000000000000, 0.774596669241483};
39 
40  /*p= 7, npoint= 4*/
41  static IssmPDouble wgt4[]={0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454};
42  static IssmPDouble xi4[]={-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053};
43 
44  /*npoint= 5*/
45  static IssmPDouble wgt5[]={0.5688888888888889,0.4786286704993665,0.4786286704993665,0.2369268850561891,0.2369268850561891};
46  static IssmPDouble xi5[]={0.0000000000000000,-0.5384693101056831,0.5384693101056831,-0.9061798459386640,0.9061798459386640};
47 
48  /*npoint= 6*/
49  static IssmPDouble wgt6[]={0.3607615730481386,0.3607615730481386,0.4679139345726910,0.4679139345726910,0.1713244923791704,0.1713244923791704};
50  static IssmPDouble xi6[]={-0.6612093864662645,0.6612093864662645,-0.2386191860831969,0.2386191860831969,-0.9324695142031521,0.9324695142031521};
51 
52  /*npoint= 7*/
53  static IssmPDouble wgt7[]={0.4179591836734694,0.3818300505051189,0.3818300505051189,0.2797053914892766,0.2797053914892766,0.1294849661688697,0.1294849661688697};
54  static IssmPDouble xi7[]={0.0000000000000000,-0.4058451513773972,0.4058451513773972,-0.7415311855993945,0.7415311855993945,-0.9491079123427585,0.9491079123427585};
55 
56  static IssmPDouble* wgtp[MAX_LINE_GAUS_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4, wgt5, wgt6, wgt7};
57  static IssmPDouble* xip [MAX_LINE_GAUS_PTS]={xi1 ,xi2 ,xi3 ,xi4, xi5, xi6, xi7};
58 
59  *pxgaus =xNew<IssmPDouble>(ngaus);
60  *pxwgt =xNew<IssmPDouble>(ngaus);
61 
62  /* check to see if Gauss points need to be calculated */
63  if(ngaus <= MAX_LINE_GAUS_PTS) {
64 
65  /* copy the points from the static arrays (noting that the pointers
66  could be returned directly, but then the calling function would
67  have to know to not free them) */
68 
69  for(i=0; i<ngaus; i++) {
70  (*pxgaus)[i]=xip [ngaus-1][i];
71  (*pxwgt )[i]=wgtp[ngaus-1][i];
72  }
73  }
74  else{
75 
76  /* calculate the Gauss points using recurrence relations */
77  alpha=xNew<IssmPDouble>(ngaus);
78  beta =xNew<IssmPDouble>(ngaus);
79 
80  /* calculate the Legendre recurrence coefficients */
81  alpha[0]=0.;
82  beta [0]=2.;
83 
84  for (i=1; i<ngaus; i++) {
85  alpha[i]=0.;
86  beta [i]=1./(4.-1./(i*i));
87  }
88 
89  /* calculate the Gauss points */
90  GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
91  xDelete<IssmPDouble>(beta);
92  xDelete<IssmPDouble>(alpha);
93  }
94 }/*}}}*/

◆ GaussLegendreTria()

void GaussLegendreTria ( int *  pngaus,
IssmDouble **  pl1,
IssmDouble **  pl2,
IssmDouble **  pl3,
IssmDouble **  pwgt,
int  iord 
)

Definition at line 95 of file GaussPoints.cpp.

95  {/*{{{*/
96  /*Gauss quadrature points for the triangle.
97 
98  Higher-order points from D.A. Dunavant, "High Degree Efficient
99  Symmetrical Gaussian Quadrature Rules for the Triangle", IJNME,
100  Vol. 21, pp. 1129-1148 (1985), as transcribed for Probe rules3.*/
101  _assert_(iord>0);
102 
103  /*Intermediaries*/
104  int i,j,ipt,nigaus;
105  IssmPDouble xi,eta;
106  IssmPDouble *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt;
107 
108  /*Hardcoded Gauss points declaration*/
109  /*p= 1, npoint= 1{{{*/
110  static IssmPDouble wgt1[]={
111  1.732050807568877};
112  static IssmPDouble l11[]={
113  0.333333333333333};
114  static IssmPDouble l21[]={
115  0.333333333333333};
116  static IssmPDouble l31[]={
117  0.333333333333333};
118  /*}}}*/
119  /*p= 2, npoint= 3 {{{*/
120  static IssmPDouble wgt2[]={
121  0.577350269189625, 0.577350269189625, 0.577350269189625};
122  static IssmPDouble l12[]={
123  0.666666666666667, 0.166666666666667, 0.166666666666667};
124  static IssmPDouble l22[]={
125  0.166666666666667, 0.666666666666667, 0.166666666666667};
126  static IssmPDouble l32[]={
127  0.166666666666667, 0.166666666666667, 0.666666666666667};
128  /*}}}*/
129  /*p= 3, npoint= 4 {{{*/
130  static IssmPDouble wgt3[]={
131  -0.974278579257493, 0.902109795608790, 0.902109795608790,
132  0.902109795608790};
133  static IssmPDouble l13[]={
134  0.333333333333333, 0.600000000000000, 0.200000000000000,
135  0.200000000000000};
136  static IssmPDouble l23[]={
137  0.333333333333333, 0.200000000000000, 0.600000000000000,
138  0.200000000000000};
139  static IssmPDouble l33[]={
140  0.333333333333333, 0.200000000000000, 0.200000000000000,
141  0.600000000000000};
142  /*}}}*/
143  /*p= 4, npoint= 6 {{{*/
144  static IssmPDouble wgt4[]={
145  0.386908262797819, 0.386908262797819, 0.386908262797819,
146  0.190442006391807, 0.190442006391807, 0.190442006391807};
147  static IssmPDouble l14[]={
148  0.108103018168070, 0.445948490915965, 0.445948490915965,
149  0.816847572980459, 0.091576213509771, 0.091576213509771};
150  static IssmPDouble l24[]={
151  0.445948490915965, 0.108103018168070, 0.445948490915965,
152  0.091576213509771, 0.816847572980459, 0.091576213509771};
153  static IssmPDouble l34[]={
154  0.445948490915965, 0.445948490915965, 0.108103018168070,
155  0.091576213509771, 0.091576213509771, 0.816847572980459};
156  /*}}}*/
157  /*p= 5, npoint= 7 {{{*/
158  static IssmPDouble wgt5[]={
159  0.389711431702997, 0.229313399254729, 0.229313399254729,
160  0.229313399254729, 0.218133059367230, 0.218133059367230,
161  0.218133059367230};
162  static IssmPDouble l15[]={
163  0.333333333333333, 0.059715871789770, 0.470142064105115,
164  0.470142064105115, 0.797426985353087, 0.101286507323456,
165  0.101286507323456};
166  static IssmPDouble l25[]={
167  0.333333333333333, 0.470142064105115, 0.059715871789770,
168  0.470142064105115, 0.101286507323456, 0.797426985353087,
169  0.101286507323456};
170  static IssmPDouble l35[]={
171  0.333333333333333, 0.470142064105115, 0.470142064105115,
172  0.059715871789770, 0.101286507323456, 0.101286507323456,
173  0.797426985353087};
174  /*}}}*/
175  /*p= 6, npoint=12 {{{*/
176  static IssmPDouble wgt6[]={
177  0.202279763184836, 0.202279763184836, 0.202279763184836,
178  0.088065961139281, 0.088065961139281, 0.088065961139281,
179  0.143502272432755, 0.143502272432755, 0.143502272432755,
180  0.143502272432755, 0.143502272432755, 0.143502272432755};
181  static IssmPDouble l16[]={
182  0.501426509658179, 0.249286745170910, 0.249286745170910,
183  0.873821971016996, 0.063089014491502, 0.063089014491502,
184  0.053145049844817, 0.053145049844817, 0.310352451033784,
185  0.636502499121399, 0.310352451033784, 0.636502499121399};
186  static IssmPDouble l26[]={
187  0.249286745170910, 0.501426509658179, 0.249286745170910,
188  0.063089014491502, 0.873821971016996, 0.063089014491502,
189  0.310352451033784, 0.636502499121399, 0.053145049844817,
190  0.053145049844817, 0.636502499121399, 0.310352451033784};
191  static IssmPDouble l36[]={
192  0.249286745170910, 0.249286745170910, 0.501426509658179,
193  0.063089014491502, 0.063089014491502, 0.873821971016996,
194  0.636502499121399, 0.310352451033784, 0.636502499121399,
195  0.310352451033784, 0.053145049844817, 0.053145049844817};
196  /*}}}*/
197  /*p= 7, npoint=13 {{{*/
198  static IssmPDouble wgt7[]={
199  -0.259062916308362, 0.304174548458604, 0.304174548458604,
200  0.304174548458604, 0.092400122517855, 0.092400122517855,
201  0.092400122517855, 0.133564951824643, 0.133564951824643,
202  0.133564951824643, 0.133564951824643, 0.133564951824643,
203  0.133564951824643};
204  static IssmPDouble l17[]={
205  0.333333333333333, 0.479308067841920, 0.260345966079040,
206  0.260345966079040, 0.869739794195568, 0.065130102902216,
207  0.065130102902216, 0.048690315425316, 0.048690315425316,
208  0.312865496004874, 0.638444188569810, 0.312865496004874,
209  0.638444188569810};
210  static IssmPDouble l27[]={
211  0.333333333333333, 0.260345966079040, 0.479308067841920,
212  0.260345966079040, 0.065130102902216, 0.869739794195568,
213  0.065130102902216, 0.312865496004874, 0.638444188569810,
214  0.048690315425316, 0.048690315425316, 0.638444188569810,
215  0.312865496004874};
216  static IssmPDouble l37[]={
217  0.333333333333333, 0.260345966079040, 0.260345966079040,
218  0.479308067841920, 0.065130102902216, 0.065130102902216,
219  0.869739794195568, 0.638444188569810, 0.312865496004874,
220  0.638444188569810, 0.312865496004874, 0.048690315425316,
221  0.048690315425316};
222  /*}}}*/
223  /*p= 8, npoint=16 {{{*/
224  static IssmPDouble wgt8[]={
225  0.249961964823104, 0.164703541925695, 0.164703541925695,
226  0.164703541925695, 0.178777729989794, 0.178777729989794,
227  0.178777729989794, 0.056219767020733, 0.056219767020733,
228  0.056219767020733, 0.047164287656184, 0.047164287656184,
229  0.047164287656184, 0.047164287656184, 0.047164287656184,
230  0.047164287656184};
231  static IssmPDouble l18[]={
232  0.333333333333333, 0.081414823414554, 0.459292588292723,
233  0.459292588292723, 0.658861384496480, 0.170569307751760,
234  0.170569307751760, 0.898905543365938, 0.050547228317031,
235  0.050547228317031, 0.008394777409958, 0.008394777409958,
236  0.263112829634638, 0.728492392955404, 0.263112829634638,
237  0.728492392955404};
238  static IssmPDouble l28[]={
239  0.333333333333333, 0.459292588292723, 0.081414823414554,
240  0.459292588292723, 0.170569307751760, 0.658861384496480,
241  0.170569307751760, 0.050547228317031, 0.898905543365938,
242  0.050547228317031, 0.263112829634638, 0.728492392955404,
243  0.008394777409958, 0.008394777409958, 0.728492392955404,
244  0.263112829634638};
245  static IssmPDouble l38[]={
246  0.333333333333333, 0.459292588292723, 0.459292588292723,
247  0.081414823414554, 0.170569307751760, 0.170569307751760,
248  0.658861384496480, 0.050547228317031, 0.050547228317031,
249  0.898905543365938, 0.728492392955404, 0.263112829634638,
250  0.728492392955404, 0.263112829634638, 0.008394777409958,
251  0.008394777409958};
252  /*}}}*/
253  /*p= 9, npoint=19 {{{*/
254  static IssmPDouble wgt9[]={
255  0.168244134395468, 0.054273292833345, 0.054273292833345,
256  0.054273292833345, 0.134801255248419, 0.134801255248419,
257  0.134801255248419, 0.137953930529909, 0.137953930529909,
258  0.137953930529909, 0.044301833780383, 0.044301833780383,
259  0.044301833780383, 0.074969289332873, 0.074969289332873,
260  0.074969289332873, 0.074969289332873, 0.074969289332873,
261  0.074969289332873};
262  static IssmPDouble l19[]={
263  0.333333333333333, 0.020634961602525, 0.489682519198738,
264  0.489682519198738, 0.125820817014127, 0.437089591492937,
265  0.437089591492937, 0.623592928761935, 0.188203535619033,
266  0.188203535619033, 0.910540973211095, 0.044729513394453,
267  0.044729513394453, 0.036838412054736, 0.036838412054736,
268  0.221962989160766, 0.741198598784498, 0.221962989160766,
269  0.741198598784498};
270  static IssmPDouble l29[]={
271  0.333333333333333, 0.489682519198738, 0.020634961602525,
272  0.489682519198738, 0.437089591492937, 0.125820817014127,
273  0.437089591492937, 0.188203535619033, 0.623592928761935,
274  0.188203535619033, 0.044729513394453, 0.910540973211095,
275  0.044729513394453, 0.221962989160766, 0.741198598784498,
276  0.036838412054736, 0.036838412054736, 0.741198598784498,
277  0.221962989160766};
278  static IssmPDouble l39[]={
279  0.333333333333333, 0.489682519198738, 0.489682519198738,
280  0.020634961602525, 0.437089591492937, 0.437089591492937,
281  0.125820817014127, 0.188203535619033, 0.188203535619033,
282  0.623592928761935, 0.044729513394453, 0.044729513394453,
283  0.910540973211095, 0.741198598784498, 0.221962989160766,
284  0.741198598784498, 0.221962989160766, 0.036838412054736,
285  0.036838412054736};
286  /*}}}*/
287  /*p=10, npoint=25 {{{*/
288  static IssmPDouble wgt10[]={
289  0.157301373584232, 0.063611224790829, 0.063611224790829,
290  0.063611224790829, 0.078498377595183, 0.078498377595183,
291  0.078498377595183, 0.126020408629139, 0.126020408629139,
292  0.126020408629139, 0.126020408629139, 0.126020408629139,
293  0.126020408629139, 0.049064223302117, 0.049064223302117,
294  0.049064223302117, 0.049064223302117, 0.049064223302117,
295  0.049064223302117, 0.016318805873179, 0.016318805873179,
296  0.016318805873179, 0.016318805873179, 0.016318805873179,
297  0.016318805873179};
298  static IssmPDouble l110[]={
299  0.333333333333333, 0.028844733232685, 0.485577633383657,
300  0.485577633383657, 0.781036849029926, 0.109481575485037,
301  0.109481575485037, 0.141707219414880, 0.141707219414880,
302  0.307939838764121, 0.550352941820999, 0.307939838764121,
303  0.550352941820999, 0.025003534762686, 0.025003534762686,
304  0.246672560639903, 0.728323904597411, 0.246672560639903,
305  0.728323904597411, 0.009540815400299, 0.009540815400299,
306  0.066803251012200, 0.923655933587500, 0.066803251012200,
307  0.923655933587500};
308  static IssmPDouble l210[]={
309  0.333333333333333, 0.485577633383657, 0.028844733232685,
310  0.485577633383657, 0.109481575485037, 0.781036849029926,
311  0.109481575485037, 0.307939838764121, 0.550352941820999,
312  0.141707219414880, 0.141707219414880, 0.550352941820999,
313  0.307939838764121, 0.246672560639903, 0.728323904597411,
314  0.025003534762686, 0.025003534762686, 0.728323904597411,
315  0.246672560639903, 0.066803251012200, 0.923655933587500,
316  0.009540815400299, 0.009540815400299, 0.923655933587500,
317  0.066803251012200};
318  static IssmPDouble l310[]={
319  0.333333333333333, 0.485577633383657, 0.485577633383657,
320  0.028844733232685, 0.109481575485037, 0.109481575485037,
321  0.781036849029926, 0.550352941820999, 0.307939838764121,
322  0.550352941820999, 0.307939838764121, 0.141707219414880,
323  0.141707219414880, 0.728323904597411, 0.246672560639903,
324  0.728323904597411, 0.246672560639903, 0.025003534762686,
325  0.025003534762686, 0.923655933587500, 0.066803251012200,
326  0.923655933587500, 0.066803251012200, 0.009540815400299,
327  0.009540815400299};
328  /*}}}*/
329  /*p=11, npoint=27 {{{*/
330  static IssmPDouble wgt11[]={
331  0.001605622060698, 0.001605622060698, 0.001605622060698,
332  0.133626914252765, 0.133626914252765, 0.133626914252765,
333  0.102750410879760, 0.102750410879760, 0.102750410879760,
334  0.062673462600454, 0.062673462600454, 0.062673462600454,
335  0.023659348114362, 0.023659348114362, 0.023659348114362,
336  0.090650537039958, 0.090650537039958, 0.090650537039958,
337  0.090650537039958, 0.090650537039958, 0.090650537039958,
338  0.035866718600836, 0.035866718600836, 0.035866718600836,
339  0.035866718600836, 0.035866718600836, 0.035866718600836};
340  static IssmPDouble l111[]={
341  -0.069222096541517, 0.534611048270758, 0.534611048270758,
342  0.202061394068290, 0.398969302965855, 0.398969302965855,
343  0.593380199137435, 0.203309900431282, 0.203309900431282,
344  0.761298175434837, 0.119350912282581, 0.119350912282581,
345  0.935270103777448, 0.032364948111276, 0.032364948111276,
346  0.050178138310495, 0.050178138310495, 0.356620648261293,
347  0.593201213428213, 0.356620648261293, 0.593201213428213,
348  0.021022016536166, 0.021022016536166, 0.171488980304042,
349  0.807489003159792, 0.171488980304042, 0.807489003159792};
350  static IssmPDouble l211[]={
351  0.534611048270758,-0.069222096541517, 0.534611048270758,
352  0.398969302965855, 0.202061394068290, 0.398969302965855,
353  0.203309900431282, 0.593380199137435, 0.203309900431282,
354  0.119350912282581, 0.761298175434837, 0.119350912282581,
355  0.032364948111276, 0.935270103777448, 0.032364948111276,
356  0.356620648261293, 0.593201213428213, 0.050178138310495,
357  0.050178138310495, 0.593201213428213, 0.356620648261293,
358  0.171488980304042, 0.807489003159792, 0.021022016536166,
359  0.021022016536166, 0.807489003159792, 0.171488980304042};
360  static IssmPDouble l311[]={
361  0.534611048270758, 0.534611048270758,-0.069222096541517,
362  0.398969302965855, 0.398969302965855, 0.202061394068290,
363  0.203309900431282, 0.203309900431282, 0.593380199137435,
364  0.119350912282581, 0.119350912282581, 0.761298175434837,
365  0.032364948111276, 0.032364948111276, 0.935270103777448,
366  0.593201213428213, 0.356620648261293, 0.593201213428213,
367  0.356620648261293, 0.050178138310495, 0.050178138310495,
368  0.807489003159792, 0.171488980304042, 0.807489003159792,
369  0.171488980304042, 0.021022016536166, 0.021022016536166};
370  /*}}}*/
371  /*p=12, npoint=33 {{{*/
372  static IssmPDouble wgt12[]={
373  0.044567514407799, 0.044567514407799, 0.044567514407799,
374  0.075677707051848, 0.075677707051848, 0.075677707051848,
375  0.108873638018933, 0.108873638018933, 0.108873638018933,
376  0.060268635501892, 0.060268635501892, 0.060268635501892,
377  0.010680277434033, 0.010680277434033, 0.010680277434033,
378  0.069925589232074, 0.069925589232074, 0.069925589232074,
379  0.069925589232074, 0.069925589232074, 0.069925589232074,
380  0.038723067079683, 0.038723067079683, 0.038723067079683,
381  0.038723067079683, 0.038723067079683, 0.038723067079683,
382  0.029992592075802, 0.029992592075802, 0.029992592075802,
383  0.029992592075802, 0.029992592075802, 0.029992592075802};
384  static IssmPDouble l112[]={
385  0.023565220452390, 0.488217389773805, 0.488217389773805,
386  0.120551215411079, 0.439724392294460, 0.439724392294460,
387  0.457579229975768, 0.271210385012116, 0.271210385012116,
388  0.744847708916828, 0.127576145541586, 0.127576145541586,
389  0.957365299093579, 0.021317350453210, 0.021317350453210,
390  0.115343494534698, 0.115343494534698, 0.275713269685514,
391  0.608943235779788, 0.275713269685514, 0.608943235779788,
392  0.022838332222257, 0.022838332222257, 0.281325580989940,
393  0.695836086787803, 0.281325580989940, 0.695836086787803,
394  0.025734050548330, 0.025734050548330, 0.116251915907597,
395  0.858014033544073, 0.116251915907597, 0.858014033544073};
396  static IssmPDouble l212[]={
397  0.488217389773805, 0.023565220452390, 0.488217389773805,
398  0.439724392294460, 0.120551215411079, 0.439724392294460,
399  0.271210385012116, 0.457579229975768, 0.271210385012116,
400  0.127576145541586, 0.744847708916828, 0.127576145541586,
401  0.021317350453210, 0.957365299093579, 0.021317350453210,
402  0.275713269685514, 0.608943235779788, 0.115343494534698,
403  0.115343494534698, 0.608943235779788, 0.275713269685514,
404  0.281325580989940, 0.695836086787803, 0.022838332222257,
405  0.022838332222257, 0.695836086787803, 0.281325580989940,
406  0.116251915907597, 0.858014033544073, 0.025734050548330,
407  0.025734050548330, 0.858014033544073, 0.116251915907597};
408  static IssmPDouble l312[]={
409  0.488217389773805, 0.488217389773805, 0.023565220452390,
410  0.439724392294460, 0.439724392294460, 0.120551215411079,
411  0.271210385012116, 0.271210385012116, 0.457579229975768,
412  0.127576145541586, 0.127576145541586, 0.744847708916828,
413  0.021317350453210, 0.021317350453210, 0.957365299093579,
414  0.608943235779788, 0.275713269685514, 0.608943235779788,
415  0.275713269685514, 0.115343494534698, 0.115343494534698,
416  0.695836086787803, 0.281325580989940, 0.695836086787803,
417  0.281325580989940, 0.022838332222257, 0.022838332222257,
418  0.858014033544073, 0.116251915907597, 0.858014033544073,
419  0.116251915907597, 0.025734050548330, 0.025734050548330};
420  /*}}}*/
421  /*p=13, npoint=37 {{{*/
422  static IssmPDouble wgt13[]={
423  0.090968907790622, 0.019537784619314, 0.019537784619314,
424  0.019537784619314, 0.054427130356344, 0.054427130356344,
425  0.054427130356344, 0.081531965976677, 0.081531965976677,
426  0.081531965976677, 0.082036138309652, 0.082036138309652,
427  0.082036138309652, 0.053983743853694, 0.053983743853694,
428  0.053983743853694, 0.013814441407066, 0.013814441407066,
429  0.013814441407066, 0.063823305703923, 0.063823305703923,
430  0.063823305703923, 0.063823305703923, 0.063823305703923,
431  0.063823305703923, 0.030140218568265, 0.030140218568265,
432  0.030140218568265, 0.030140218568265, 0.030140218568265,
433  0.030140218568265, 0.026884523429480, 0.026884523429480,
434  0.026884523429480, 0.026884523429480, 0.026884523429480,
435  0.026884523429480};
436  static IssmPDouble l113[]={
437  0.333333333333333, 0.009903630120591, 0.495048184939705,
438  0.495048184939705, 0.062566729780852, 0.468716635109574,
439  0.468716635109574, 0.170957326397447, 0.414521336801277,
440  0.414521336801277, 0.541200855914337, 0.229399572042831,
441  0.229399572042831, 0.771151009607340, 0.114424495196330,
442  0.114424495196330, 0.950377217273082, 0.024811391363459,
443  0.024811391363459, 0.094853828379579, 0.094853828379579,
444  0.268794997058761, 0.636351174561660, 0.268794997058761,
445  0.636351174561660, 0.018100773278807, 0.018100773278807,
446  0.291730066734288, 0.690169159986905, 0.291730066734288,
447  0.690169159986905, 0.022233076674090, 0.022233076674090,
448  0.126357385491669, 0.851409537834241, 0.126357385491669,
449  0.851409537834241};
450  static IssmPDouble l213[]={
451  0.333333333333333, 0.495048184939705, 0.009903630120591,
452  0.495048184939705, 0.468716635109574, 0.062566729780852,
453  0.468716635109574, 0.414521336801277, 0.170957326397447,
454  0.414521336801277, 0.229399572042831, 0.541200855914337,
455  0.229399572042831, 0.114424495196330, 0.771151009607340,
456  0.114424495196330, 0.024811391363459, 0.950377217273082,
457  0.024811391363459, 0.268794997058761, 0.636351174561660,
458  0.094853828379579, 0.094853828379579, 0.636351174561660,
459  0.268794997058761, 0.291730066734288, 0.690169159986905,
460  0.018100773278807, 0.018100773278807, 0.690169159986905,
461  0.291730066734288, 0.126357385491669, 0.851409537834241,
462  0.022233076674090, 0.022233076674090, 0.851409537834241,
463  0.126357385491669};
464  static IssmPDouble l313[]={
465  0.333333333333333, 0.495048184939705, 0.495048184939705,
466  0.009903630120591, 0.468716635109574, 0.468716635109574,
467  0.062566729780852, 0.414521336801277, 0.414521336801277,
468  0.170957326397447, 0.229399572042831, 0.229399572042831,
469  0.541200855914337, 0.114424495196330, 0.114424495196330,
470  0.771151009607340, 0.024811391363459, 0.024811391363459,
471  0.950377217273082, 0.636351174561660, 0.268794997058761,
472  0.636351174561660, 0.268794997058761, 0.094853828379579,
473  0.094853828379579, 0.690169159986905, 0.291730066734288,
474  0.690169159986905, 0.291730066734288, 0.018100773278807,
475  0.018100773278807, 0.851409537834241, 0.126357385491669,
476  0.851409537834241, 0.126357385491669, 0.022233076674090,
477  0.022233076674090};
478  /*}}}*/
479  /*p=14, npoint=42{{{*/
480  static IssmPDouble wgt14[]={
481  0.037903474783419, 0.037903474783419, 0.037903474783419,
482  0.056791094234956, 0.056791094234956, 0.056791094234956,
483  0.089675379523011, 0.089675379523011, 0.089675379523011,
484  0.073027745871103, 0.073027745871103, 0.073027745871103,
485  0.024999901169244, 0.024999901169244, 0.024999901169244,
486  0.008527585185524, 0.008527585185524, 0.008527585185524,
487  0.042722337771116, 0.042722337771116, 0.042722337771116,
488  0.042722337771116, 0.042722337771116, 0.042722337771116,
489  0.066807816407881, 0.066807816407881, 0.066807816407881,
490  0.066807816407881, 0.066807816407881, 0.066807816407881,
491  0.025004419126360, 0.025004419126360, 0.025004419126360,
492  0.025004419126360, 0.025004419126360, 0.025004419126360,
493  0.008677970905831, 0.008677970905831, 0.008677970905831,
494  0.008677970905831, 0.008677970905831, 0.008677970905831};
495  static IssmPDouble l114[]={
496  0.022072179275643, 0.488963910362179, 0.488963910362179,
497  0.164710561319092, 0.417644719340454, 0.417644719340454,
498  0.453044943382323, 0.273477528308839, 0.273477528308839,
499  0.645588935174913, 0.177205532412543, 0.177205532412543,
500  0.876400233818255, 0.061799883090873, 0.061799883090873,
501  0.961218077502598, 0.019390961248701, 0.019390961248701,
502  0.057124757403648, 0.057124757403648, 0.172266687821356,
503  0.770608554774996, 0.172266687821356, 0.770608554774996,
504  0.092916249356972, 0.092916249356972, 0.336861459796345,
505  0.570222290846683, 0.336861459796345, 0.570222290846683,
506  0.014646950055654, 0.014646950055654, 0.298372882136258,
507  0.686980167808088, 0.298372882136258, 0.686980167808088,
508  0.001268330932872, 0.001268330932872, 0.118974497696957,
509  0.879757171370171, 0.118974497696957, 0.879757171370171};
510  static IssmPDouble l214[]={
511  0.488963910362179, 0.022072179275643, 0.488963910362179,
512  0.417644719340454, 0.164710561319092, 0.417644719340454,
513  0.273477528308839, 0.453044943382323, 0.273477528308839,
514  0.177205532412543, 0.645588935174913, 0.177205532412543,
515  0.061799883090873, 0.876400233818255, 0.061799883090873,
516  0.019390961248701, 0.961218077502598, 0.019390961248701,
517  0.172266687821356, 0.770608554774996, 0.057124757403648,
518  0.057124757403648, 0.770608554774996, 0.172266687821356,
519  0.336861459796345, 0.570222290846683, 0.092916249356972,
520  0.092916249356972, 0.570222290846683, 0.336861459796345,
521  0.298372882136258, 0.686980167808088, 0.014646950055654,
522  0.014646950055654, 0.686980167808088, 0.298372882136258,
523  0.118974497696957, 0.879757171370171, 0.001268330932872,
524  0.001268330932872, 0.879757171370171, 0.118974497696957};
525  static IssmPDouble l314[]={
526  0.488963910362179, 0.488963910362179, 0.022072179275643,
527  0.417644719340454, 0.417644719340454, 0.164710561319092,
528  0.273477528308839, 0.273477528308839, 0.453044943382323,
529  0.177205532412543, 0.177205532412543, 0.645588935174913,
530  0.061799883090873, 0.061799883090873, 0.876400233818255,
531  0.019390961248701, 0.019390961248701, 0.961218077502598,
532  0.770608554774996, 0.172266687821356, 0.770608554774996,
533  0.172266687821356, 0.057124757403648, 0.057124757403648,
534  0.570222290846683, 0.336861459796345, 0.570222290846683,
535  0.336861459796345, 0.092916249356972, 0.092916249356972,
536  0.686980167808088, 0.298372882136258, 0.686980167808088,
537  0.298372882136258, 0.014646950055654, 0.014646950055654,
538  0.879757171370171, 0.118974497696957, 0.879757171370171,
539  0.118974497696957, 0.001268330932872, 0.001268330932872};
540  /*}}}*/
541  /*p=15, npoint=48{{{*/
542  static IssmPDouble wgt15[]={
543  0.003320126005206, 0.003320126005206, 0.003320126005206,
544  0.076641563419124, 0.076641563419124, 0.076641563419124,
545  0.088657703045151, 0.088657703045151, 0.088657703045151,
546  0.041028362044303, 0.041028362044303, 0.041028362044303,
547  0.023018566716310, 0.023018566716310, 0.023018566716310,
548  0.008225364846296, 0.008225364846296, 0.008225364846296,
549  0.066770684377964, 0.066770684377964, 0.066770684377964,
550  0.066770684377964, 0.066770684377964, 0.066770684377964,
551  0.047139173172681, 0.047139173172681, 0.047139173172681,
552  0.047139173172681, 0.047139173172681, 0.047139173172681,
553  0.003779468865339, 0.003779468865339, 0.003779468865339,
554  0.003779468865339, 0.003779468865339, 0.003779468865339,
555  0.037248306609289, 0.037248306609289, 0.037248306609289,
556  0.037248306609289, 0.037248306609289, 0.037248306609289,
557  0.013291658531346, 0.013291658531346, 0.013291658531346,
558  0.013291658531346, 0.013291658531346, 0.013291658531346};
559  static IssmPDouble l115[]={
560  -0.013945833716486, 0.506972916858243, 0.506972916858243,
561  0.137187291433955, 0.431406354283023, 0.431406354283023,
562  0.444612710305711, 0.277693644847144, 0.277693644847144,
563  0.747070217917492, 0.126464891041254, 0.126464891041254,
564  0.858383228050628, 0.070808385974686, 0.070808385974686,
565  0.962069659517853, 0.018965170241073, 0.018965170241073,
566  0.133734161966621, 0.133734161966621, 0.261311371140087,
567  0.604954466893291, 0.261311371140087, 0.604954466893291,
568  0.036366677396917, 0.036366677396917, 0.388046767090269,
569  0.575586555512814, 0.388046767090269, 0.575586555512814,
570  -0.010174883126571,-0.010174883126571, 0.285712220049916,
571  0.724462663076655, 0.285712220049916, 0.724462663076655,
572  0.036843869875878, 0.036843869875878, 0.215599664072284,
573  0.747556466051838, 0.215599664072284, 0.747556466051838,
574  0.012459809331199, 0.012459809331199, 0.103575616576386,
575  0.883964574092416, 0.103575616576386, 0.883964574092416};
576  static IssmPDouble l215[]={
577  0.506972916858243,-0.013945833716486, 0.506972916858243,
578  0.431406354283023, 0.137187291433955, 0.431406354283023,
579  0.277693644847144, 0.444612710305711, 0.277693644847144,
580  0.126464891041254, 0.747070217917492, 0.126464891041254,
581  0.070808385974686, 0.858383228050628, 0.070808385974686,
582  0.018965170241073, 0.962069659517853, 0.018965170241073,
583  0.261311371140087, 0.604954466893291, 0.133734161966621,
584  0.133734161966621, 0.604954466893291, 0.261311371140087,
585  0.388046767090269, 0.575586555512814, 0.036366677396917,
586  0.036366677396917, 0.575586555512814, 0.388046767090269,
587  0.285712220049916, 0.724462663076655,-0.010174883126571,
588  -0.010174883126571, 0.724462663076655, 0.285712220049916,
589  0.215599664072284, 0.747556466051838, 0.036843869875878,
590  0.036843869875878, 0.747556466051838, 0.215599664072284,
591  0.103575616576386, 0.883964574092416, 0.012459809331199,
592  0.012459809331199, 0.883964574092416, 0.103575616576386};
593  static IssmPDouble l315[]={
594  0.506972916858243, 0.506972916858243,-0.013945833716486,
595  0.431406354283023, 0.431406354283023, 0.137187291433955,
596  0.277693644847144, 0.277693644847144, 0.444612710305711,
597  0.126464891041254, 0.126464891041254, 0.747070217917492,
598  0.070808385974686, 0.070808385974686, 0.858383228050628,
599  0.018965170241073, 0.018965170241073, 0.962069659517853,
600  0.604954466893291, 0.261311371140087, 0.604954466893291,
601  0.261311371140087, 0.133734161966621, 0.133734161966621,
602  0.575586555512814, 0.388046767090269, 0.575586555512814,
603  0.388046767090269, 0.036366677396917, 0.036366677396917,
604  0.724462663076655, 0.285712220049916, 0.724462663076655,
605  0.285712220049916,-0.010174883126571,-0.010174883126571,
606  0.747556466051838, 0.215599664072284, 0.747556466051838,
607  0.215599664072284, 0.036843869875878, 0.036843869875878,
608  0.883964574092416, 0.103575616576386, 0.883964574092416,
609  0.103575616576386, 0.012459809331199, 0.012459809331199};
610  /*}}}*/
611  /*p=16, npoint=52 {{{*/
612  static IssmPDouble wgt16[]={
613  0.081191089584902, 0.011095307165226, 0.011095307165226,
614  0.011095307165226, 0.072244353151393, 0.072244353151393,
615  0.072244353151393, 0.046577417012049, 0.046577417012049,
616  0.046577417012049, 0.072975670074230, 0.072975670074230,
617  0.072975670074230, 0.051961986412307, 0.051961986412307,
618  0.051961986412307, 0.024595292810646, 0.024595292810646,
619  0.024595292810646, 0.006205006808607, 0.006205006808607,
620  0.006205006808607, 0.056764756525753, 0.056764756525753,
621  0.056764756525753, 0.056764756525753, 0.056764756525753,
622  0.056764756525753, 0.026497443692048, 0.026497443692048,
623  0.026497443692048, 0.026497443692048, 0.026497443692048,
624  0.026497443692048, 0.004133096181263, 0.004133096181263,
625  0.004133096181263, 0.004133096181263, 0.004133096181263,
626  0.004133096181263, 0.033055830705140, 0.033055830705140,
627  0.033055830705140, 0.033055830705140, 0.033055830705140,
628  0.033055830705140, 0.011864642509229, 0.011864642509229,
629  0.011864642509229, 0.011864642509229, 0.011864642509229,
630  0.011864642509229};
631  static IssmPDouble l116[]={
632  0.333333333333333, 0.005238916103123, 0.497380541948438,
633  0.497380541948438, 0.173061122901295, 0.413469438549352,
634  0.413469438549352, 0.059082801866017, 0.470458599066991,
635  0.470458599066991, 0.518892500060958, 0.240553749969521,
636  0.240553749969521, 0.704068411554854, 0.147965794222573,
637  0.147965794222573, 0.849069624685052, 0.075465187657474,
638  0.075465187657474, 0.966807194753950, 0.016596402623025,
639  0.016596402623025, 0.103575692245252, 0.103575692245252,
640  0.296555596579887, 0.599868711174861, 0.296555596579887,
641  0.599868711174861, 0.020083411655416, 0.020083411655416,
642  0.337723063403079, 0.642193524941505, 0.337723063403079,
643  0.642193524941505,-0.004341002614139,-0.004341002614139,
644  0.204748281642812, 0.799592720971327, 0.204748281642812,
645  0.799592720971327, 0.041941786468010, 0.041941786468010,
646  0.189358492130623, 0.768699721401368, 0.189358492130623,
647  0.768699721401368, 0.014317320230681, 0.014317320230681,
648  0.085283615682657, 0.900399064086661, 0.085283615682657,
649  0.900399064086661};
650  static IssmPDouble l216[]={
651  0.333333333333333, 0.497380541948438, 0.005238916103123,
652  0.497380541948438, 0.413469438549352, 0.173061122901295,
653  0.413469438549352, 0.470458599066991, 0.059082801866017,
654  0.470458599066991, 0.240553749969521, 0.518892500060958,
655  0.240553749969521, 0.147965794222573, 0.704068411554854,
656  0.147965794222573, 0.075465187657474, 0.849069624685052,
657  0.075465187657474, 0.016596402623025, 0.966807194753950,
658  0.016596402623025, 0.296555596579887, 0.599868711174861,
659  0.103575692245252, 0.103575692245252, 0.599868711174861,
660  0.296555596579887, 0.337723063403079, 0.642193524941505,
661  0.020083411655416, 0.020083411655416, 0.642193524941505,
662  0.337723063403079, 0.204748281642812, 0.799592720971327,
663  -0.004341002614139,-0.004341002614139, 0.799592720971327,
664  0.204748281642812, 0.189358492130623, 0.768699721401368,
665  0.041941786468010, 0.041941786468010, 0.768699721401368,
666  0.189358492130623, 0.085283615682657, 0.900399064086661,
667  0.014317320230681, 0.014317320230681, 0.900399064086661,
668  0.085283615682657};
669  static IssmPDouble l316[]={
670  0.333333333333333, 0.497380541948438, 0.497380541948438,
671  0.005238916103123, 0.413469438549352, 0.413469438549352,
672  0.173061122901295, 0.470458599066991, 0.470458599066991,
673  0.059082801866017, 0.240553749969521, 0.240553749969521,
674  0.518892500060958, 0.147965794222573, 0.147965794222573,
675  0.704068411554854, 0.075465187657474, 0.075465187657474,
676  0.849069624685052, 0.016596402623025, 0.016596402623025,
677  0.966807194753950, 0.599868711174861, 0.296555596579887,
678  0.599868711174861, 0.296555596579887, 0.103575692245252,
679  0.103575692245252, 0.642193524941505, 0.337723063403079,
680  0.642193524941505, 0.337723063403079, 0.020083411655416,
681  0.020083411655416, 0.799592720971327, 0.204748281642812,
682  0.799592720971327, 0.204748281642812,-0.004341002614139,
683  -0.004341002614139, 0.768699721401368, 0.189358492130623,
684  0.768699721401368, 0.189358492130623, 0.041941786468010,
685  0.041941786468010, 0.900399064086661, 0.085283615682657,
686  0.900399064086661, 0.085283615682657, 0.014317320230681,
687  0.014317320230681};
688  /*}}}*/
689  /*p=17, npoint=61{{{*/
690  static IssmPDouble wgt17[]={
691  0.057914928034477, 0.008822054327014, 0.008822054327014,
692  0.008822054327014, 0.025410682752829, 0.025410682752829,
693  0.025410682752829, 0.042176958517489, 0.042176958517489,
694  0.042176958517489, 0.053879858604088, 0.053879858604088,
695  0.053879858604088, 0.054138904728481, 0.054138904728481,
696  0.054138904728481, 0.042981974139367, 0.042981974139367,
697  0.042981974139367, 0.024345832713105, 0.024345832713105,
698  0.024345832713105, 0.005533341446715, 0.005533341446715,
699  0.005533341446715, 0.014063655552443, 0.014063655552443,
700  0.014063655552443, 0.014063655552443, 0.014063655552443,
701  0.014063655552443, 0.046428907569036, 0.046428907569036,
702  0.046428907569036, 0.046428907569036, 0.046428907569036,
703  0.046428907569036, 0.031973646148520, 0.031973646148520,
704  0.031973646148520, 0.031973646148520, 0.031973646148520,
705  0.031973646148520, 0.014682366990538, 0.014682366990538,
706  0.014682366990538, 0.014682366990538, 0.014682366990538,
707  0.014682366990538, 0.031684053418215, 0.031684053418215,
708  0.031684053418215, 0.031684053418215, 0.031684053418215,
709  0.031684053418215, 0.011545213295771, 0.011545213295771,
710  0.011545213295771, 0.011545213295771, 0.011545213295771,
711  0.011545213295771};
712  static IssmPDouble l117[]={
713  0.333333333333333, 0.005658918886452, 0.497170540556774,
714  0.497170540556774, 0.035647354750751, 0.482176322624625,
715  0.482176322624625, 0.099520061958437, 0.450239969020782,
716  0.450239969020782, 0.199467521245206, 0.400266239377397,
717  0.400266239377397, 0.495717464058095, 0.252141267970953,
718  0.252141267970953, 0.675905990683077, 0.162047004658461,
719  0.162047004658461, 0.848248235478508, 0.075875882260746,
720  0.075875882260746, 0.968690546064356, 0.015654726967822,
721  0.015654726967822, 0.010186928826919, 0.010186928826919,
722  0.334319867363658, 0.655493203809423, 0.334319867363658,
723  0.655493203809423, 0.135440871671036, 0.135440871671036,
724  0.292221537796944, 0.572337590532020, 0.292221537796944,
725  0.572337590532020, 0.054423924290583, 0.054423924290583,
726  0.319574885423190, 0.626001190286228, 0.319574885423190,
727  0.626001190286228, 0.012868560833637, 0.012868560833637,
728  0.190704224192292, 0.796427214974071, 0.190704224192292,
729  0.796427214974071, 0.067165782413524, 0.067165782413524,
730  0.180483211648746, 0.752351005937729, 0.180483211648746,
731  0.752351005937729, 0.014663182224828, 0.014663182224828,
732  0.080711313679564, 0.904625504095608, 0.080711313679564,
733  0.904625504095608};
734  static IssmPDouble l217[]={
735  0.333333333333333, 0.497170540556774, 0.005658918886452,
736  0.497170540556774, 0.482176322624625, 0.035647354750751,
737  0.482176322624625, 0.450239969020782, 0.099520061958437,
738  0.450239969020782, 0.400266239377397, 0.199467521245206,
739  0.400266239377397, 0.252141267970953, 0.495717464058095,
740  0.252141267970953, 0.162047004658461, 0.675905990683077,
741  0.162047004658461, 0.075875882260746, 0.848248235478508,
742  0.075875882260746, 0.015654726967822, 0.968690546064356,
743  0.015654726967822, 0.334319867363658, 0.655493203809423,
744  0.010186928826919, 0.010186928826919, 0.655493203809423,
745  0.334319867363658, 0.292221537796944, 0.572337590532020,
746  0.135440871671036, 0.135440871671036, 0.572337590532020,
747  0.292221537796944, 0.319574885423190, 0.626001190286228,
748  0.054423924290583, 0.054423924290583, 0.626001190286228,
749  0.319574885423190, 0.190704224192292, 0.796427214974071,
750  0.012868560833637, 0.012868560833637, 0.796427214974071,
751  0.190704224192292, 0.180483211648746, 0.752351005937729,
752  0.067165782413524, 0.067165782413524, 0.752351005937729,
753  0.180483211648746, 0.080711313679564, 0.904625504095608,
754  0.014663182224828, 0.014663182224828, 0.904625504095608,
755  0.080711313679564};
756  static IssmPDouble l317[]={
757  0.333333333333333, 0.497170540556774, 0.497170540556774,
758  0.005658918886452, 0.482176322624625, 0.482176322624625,
759  0.035647354750751, 0.450239969020782, 0.450239969020782,
760  0.099520061958437, 0.400266239377397, 0.400266239377397,
761  0.199467521245206, 0.252141267970953, 0.252141267970953,
762  0.495717464058095, 0.162047004658461, 0.162047004658461,
763  0.675905990683077, 0.075875882260746, 0.075875882260746,
764  0.848248235478508, 0.015654726967822, 0.015654726967822,
765  0.968690546064356, 0.655493203809423, 0.334319867363658,
766  0.655493203809423, 0.334319867363658, 0.010186928826919,
767  0.010186928826919, 0.572337590532020, 0.292221537796944,
768  0.572337590532020, 0.292221537796944, 0.135440871671036,
769  0.135440871671036, 0.626001190286228, 0.319574885423190,
770  0.626001190286228, 0.319574885423190, 0.054423924290583,
771  0.054423924290583, 0.796427214974071, 0.190704224192292,
772  0.796427214974071, 0.190704224192292, 0.012868560833637,
773  0.012868560833637, 0.752351005937729, 0.180483211648746,
774  0.752351005937729, 0.180483211648746, 0.067165782413524,
775  0.067165782413524, 0.904625504095608, 0.080711313679564,
776  0.904625504095608, 0.080711313679564, 0.014663182224828,
777  0.014663182224828};
778  /*}}}*/
779  /*p=18, npoint=70 {{{*/
780 
781  static IssmPDouble wgt18[]={
782  0.053364381350150, 0.015713921277179, 0.015713921277179,
783  0.015713921277179, 0.032495554156279, 0.032495554156279,
784  0.032495554156279, 0.033672969465771, 0.033672969465771,
785  0.033672969465771, 0.048071249104579, 0.048071249104579,
786  0.048071249104579, 0.055869421169115, 0.055869421169115,
787  0.055869421169115, 0.043429498443148, 0.043429498443148,
788  0.043429498443148, 0.026451755176745, 0.026451755176745,
789  0.026451755176745, 0.011767418126433, 0.011767418126433,
790  0.011767418126433,-0.003850519950463,-0.003850519950463,
791  -0.003850519950463, 0.010967196889496, 0.010967196889496,
792  0.010967196889496, 0.010967196889496, 0.010967196889496,
793  0.010967196889496, 0.047211440790349, 0.047211440790349,
794  0.047211440790349, 0.047211440790349, 0.047211440790349,
795  0.047211440790349, 0.030617090859378, 0.030617090859378,
796  0.030617090859378, 0.030617090859378, 0.030617090859378,
797  0.030617090859378, 0.031834201210069, 0.031834201210069,
798  0.031834201210069, 0.031834201210069, 0.031834201210069,
799  0.031834201210069, 0.014037809005559, 0.014037809005559,
800  0.014037809005559, 0.014037809005559, 0.014037809005559,
801  0.014037809005559, 0.013222699422034, 0.013222699422034,
802  0.013222699422034, 0.013222699422034, 0.013222699422034,
803  0.013222699422034, 0.000079999375178, 0.000079999375178,
804  0.000079999375178, 0.000079999375178, 0.000079999375178,
805  0.000079999375178};
806  static IssmPDouble l118[]={
807  0.333333333333333, 0.013310382738157, 0.493344808630921,
808  0.493344808630921, 0.061578811516086, 0.469210594241957,
809  0.469210594241957, 0.127437208225989, 0.436281395887006,
810  0.436281395887006, 0.210307658653168, 0.394846170673416,
811  0.394846170673416, 0.500410862393686, 0.249794568803157,
812  0.249794568803157, 0.677135612512315, 0.161432193743843,
813  0.161432193743843, 0.846803545029257, 0.076598227485371,
814  0.076598227485371, 0.951495121293100, 0.024252439353450,
815  0.024252439353450, 0.913707265566071, 0.043146367216965,
816  0.043146367216965, 0.008430536202420, 0.008430536202420,
817  0.358911494940944, 0.632657968856636, 0.358911494940944,
818  0.632657968856636, 0.131186551737188, 0.131186551737188,
819  0.294402476751957, 0.574410971510855, 0.294402476751957,
820  0.574410971510855, 0.050203151565675, 0.050203151565675,
821  0.325017801641814, 0.624779046792512, 0.325017801641814,
822  0.624779046792512, 0.066329263810916, 0.066329263810916,
823  0.184737559666046, 0.748933176523037, 0.184737559666046,
824  0.748933176523037, 0.011996194566236, 0.011996194566236,
825  0.218796800013321, 0.769207005420443, 0.218796800013321,
826  0.769207005420443, 0.014858100590125, 0.014858100590125,
827  0.101179597136408, 0.883962302273467, 0.101179597136408,
828  0.883962302273467,-0.035222015287949,-0.035222015287949,
829  0.020874755282586, 1.014347260005363, 0.020874755282586,
830  1.014347260005363};
831  static IssmPDouble l218[]={
832  0.333333333333333, 0.493344808630921, 0.013310382738157,
833  0.493344808630921, 0.469210594241957, 0.061578811516086,
834  0.469210594241957, 0.436281395887006, 0.127437208225989,
835  0.436281395887006, 0.394846170673416, 0.210307658653168,
836  0.394846170673416, 0.249794568803157, 0.500410862393686,
837  0.249794568803157, 0.161432193743843, 0.677135612512315,
838  0.161432193743843, 0.076598227485371, 0.846803545029257,
839  0.076598227485371, 0.024252439353450, 0.951495121293100,
840  0.024252439353450, 0.043146367216965, 0.913707265566071,
841  0.043146367216965, 0.358911494940944, 0.632657968856636,
842  0.008430536202420, 0.008430536202420, 0.632657968856636,
843  0.358911494940944, 0.294402476751957, 0.574410971510855,
844  0.131186551737188, 0.131186551737188, 0.574410971510855,
845  0.294402476751957, 0.325017801641814, 0.624779046792512,
846  0.050203151565675, 0.050203151565675, 0.624779046792512,
847  0.325017801641814, 0.184737559666046, 0.748933176523037,
848  0.066329263810916, 0.066329263810916, 0.748933176523037,
849  0.184737559666046, 0.218796800013321, 0.769207005420443,
850  0.011996194566236, 0.011996194566236, 0.769207005420443,
851  0.218796800013321, 0.101179597136408, 0.883962302273467,
852  0.014858100590125, 0.014858100590125, 0.883962302273467,
853  0.101179597136408, 0.020874755282586, 1.014347260005363,
854  -0.035222015287949,-0.035222015287949, 1.014347260005363,
855  0.020874755282586};
856  static IssmPDouble l318[]={
857  0.333333333333333, 0.493344808630921, 0.493344808630921,
858  0.013310382738157, 0.469210594241957, 0.469210594241957,
859  0.061578811516086, 0.436281395887006, 0.436281395887006,
860  0.127437208225989, 0.394846170673416, 0.394846170673416,
861  0.210307658653168, 0.249794568803157, 0.249794568803157,
862  0.500410862393686, 0.161432193743843, 0.161432193743843,
863  0.677135612512315, 0.076598227485371, 0.076598227485371,
864  0.846803545029257, 0.024252439353450, 0.024252439353450,
865  0.951495121293100, 0.043146367216965, 0.043146367216965,
866  0.913707265566071, 0.632657968856636, 0.358911494940944,
867  0.632657968856636, 0.358911494940944, 0.008430536202420,
868  0.008430536202420, 0.574410971510855, 0.294402476751957,
869  0.574410971510855, 0.294402476751957, 0.131186551737188,
870  0.131186551737188, 0.624779046792512, 0.325017801641814,
871  0.624779046792512, 0.325017801641814, 0.050203151565675,
872  0.050203151565675, 0.748933176523037, 0.184737559666046,
873  0.748933176523037, 0.184737559666046, 0.066329263810916,
874  0.066329263810916, 0.769207005420443, 0.218796800013321,
875  0.769207005420443, 0.218796800013321, 0.011996194566236,
876  0.011996194566236, 0.883962302273467, 0.101179597136408,
877  0.883962302273467, 0.101179597136408, 0.014858100590125,
878  0.014858100590125, 1.014347260005363, 0.020874755282586,
879  1.014347260005363, 0.020874755282586,-0.035222015287949,
880  -0.035222015287949};
881  /*}}}*/
882  /*p=19, npoint=73 {{{*/
883 
884  static IssmPDouble wgt19[]={
885  0.056995437856306, 0.017893352515055, 0.017893352515055,
886  0.017893352515055, 0.038775849701151, 0.038775849701151,
887  0.038775849701151, 0.052422467754193, 0.052422467754193,
888  0.052422467754193, 0.052811905405354, 0.052811905405354,
889  0.052811905405354, 0.041844983939388, 0.041844983939388,
890  0.041844983939388, 0.027800807314648, 0.027800807314648,
891  0.027800807314648, 0.014002903771278, 0.014002903771278,
892  0.014002903771278, 0.003601560678933, 0.003601560678933,
893  0.003601560678933, 0.006728804180578, 0.006728804180578,
894  0.006728804180578, 0.006728804180578, 0.006728804180578,
895  0.006728804180578, 0.044295745540949, 0.044295745540949,
896  0.044295745540949, 0.044295745540949, 0.044295745540949,
897  0.044295745540949, 0.015382176206141, 0.015382176206141,
898  0.015382176206141, 0.015382176206141, 0.015382176206141,
899  0.015382176206141, 0.027928534240338, 0.027928534240338,
900  0.027928534240338, 0.027928534240338, 0.027928534240338,
901  0.027928534240338, 0.004316169837400, 0.004316169837400,
902  0.004316169837400, 0.004316169837400, 0.004316169837400,
903  0.004316169837400, 0.031597525960379, 0.031597525960379,
904  0.031597525960379, 0.031597525960379, 0.031597525960379,
905  0.031597525960379, 0.017768353603780, 0.017768353603780,
906  0.017768353603780, 0.017768353603780, 0.017768353603780,
907  0.017768353603780, 0.006581669842530, 0.006581669842530,
908  0.006581669842530, 0.006581669842530, 0.006581669842530,
909  0.006581669842530};
910  static IssmPDouble l119[]={
911  0.333333333333333, 0.020780025853987, 0.489609987073006,
912  0.489609987073006, 0.090926214604215, 0.454536892697893,
913  0.454536892697893, 0.197166638701138, 0.401416680649431,
914  0.401416680649431, 0.488896691193805, 0.255551654403098,
915  0.255551654403098, 0.645844115695741, 0.177077942152130,
916  0.177077942152130, 0.779877893544096, 0.110061053227952,
917  0.110061053227952, 0.888942751496321, 0.055528624251840,
918  0.055528624251840, 0.974756272445543, 0.012621863777229,
919  0.012621863777229, 0.003611417848412, 0.003611417848412,
920  0.395754787356943, 0.600633794794645, 0.395754787356943,
921  0.600633794794645, 0.134466754530780, 0.134466754530780,
922  0.307929983880436, 0.557603261588784, 0.307929983880436,
923  0.557603261588784, 0.014446025776115, 0.014446025776115,
924  0.264566948406520, 0.720987025817365, 0.264566948406520,
925  0.720987025817365, 0.046933578838178, 0.046933578838178,
926  0.358539352205951, 0.594527068955871, 0.358539352205951,
927  0.594527068955871, 0.002861120350567, 0.002861120350567,
928  0.157807405968595, 0.839331473680839, 0.157807405968595,
929  0.839331473680839, 0.223861424097916, 0.223861424097916,
930  0.075050596975911, 0.701087978926173, 0.075050596975911,
931  0.701087978926173, 0.034647074816760, 0.034647074816760,
932  0.142421601113383, 0.822931324069857, 0.142421601113383,
933  0.822931324069857, 0.010161119296278, 0.010161119296278,
934  0.065494628082938, 0.924344252620784, 0.065494628082938,
935  0.924344252620784};
936  static IssmPDouble l219[]={
937  0.333333333333333, 0.489609987073006, 0.020780025853987,
938  0.489609987073006, 0.454536892697893, 0.090926214604215,
939  0.454536892697893, 0.401416680649431, 0.197166638701138,
940  0.401416680649431, 0.255551654403098, 0.488896691193805,
941  0.255551654403098, 0.177077942152130, 0.645844115695741,
942  0.177077942152130, 0.110061053227952, 0.779877893544096,
943  0.110061053227952, 0.055528624251840, 0.888942751496321,
944  0.055528624251840, 0.012621863777229, 0.974756272445543,
945  0.012621863777229, 0.395754787356943, 0.600633794794645,
946  0.003611417848412, 0.003611417848412, 0.600633794794645,
947  0.395754787356943, 0.307929983880436, 0.557603261588784,
948  0.134466754530780, 0.134466754530780, 0.557603261588784,
949  0.307929983880436, 0.264566948406520, 0.720987025817365,
950  0.014446025776115, 0.014446025776115, 0.720987025817365,
951  0.264566948406520, 0.358539352205951, 0.594527068955871,
952  0.046933578838178, 0.046933578838178, 0.594527068955871,
953  0.358539352205951, 0.157807405968595, 0.839331473680839,
954  0.002861120350567, 0.002861120350567, 0.839331473680839,
955  0.157807405968595, 0.075050596975911, 0.701087978926173,
956  0.223861424097916, 0.223861424097916, 0.701087978926173,
957  0.075050596975911, 0.142421601113383, 0.822931324069857,
958  0.034647074816760, 0.034647074816760, 0.822931324069857,
959  0.142421601113383, 0.065494628082938, 0.924344252620784,
960  0.010161119296278, 0.010161119296278, 0.924344252620784,
961  0.065494628082938};
962  static IssmPDouble l319[]={
963  0.333333333333333, 0.489609987073006, 0.489609987073006,
964  0.020780025853987, 0.454536892697893, 0.454536892697893,
965  0.090926214604215, 0.401416680649431, 0.401416680649431,
966  0.197166638701138, 0.255551654403098, 0.255551654403098,
967  0.488896691193805, 0.177077942152130, 0.177077942152130,
968  0.645844115695741, 0.110061053227952, 0.110061053227952,
969  0.779877893544096, 0.055528624251840, 0.055528624251840,
970  0.888942751496321, 0.012621863777229, 0.012621863777229,
971  0.974756272445543, 0.600633794794645, 0.395754787356943,
972  0.600633794794645, 0.395754787356943, 0.003611417848412,
973  0.003611417848412, 0.557603261588784, 0.307929983880436,
974  0.557603261588784, 0.307929983880436, 0.134466754530780,
975  0.134466754530780, 0.720987025817365, 0.264566948406520,
976  0.720987025817365, 0.264566948406520, 0.014446025776115,
977  0.014446025776115, 0.594527068955871, 0.358539352205951,
978  0.594527068955871, 0.358539352205951, 0.046933578838178,
979  0.046933578838178, 0.839331473680839, 0.157807405968595,
980  0.839331473680839, 0.157807405968595, 0.002861120350567,
981  0.002861120350567, 0.701087978926173, 0.075050596975911,
982  0.701087978926173, 0.075050596975911, 0.223861424097916,
983  0.223861424097916, 0.822931324069857, 0.142421601113383,
984  0.822931324069857, 0.142421601113383, 0.034647074816760,
985  0.034647074816760, 0.924344252620784, 0.065494628082938,
986  0.924344252620784, 0.065494628082938, 0.010161119296278,
987  0.010161119296278};
988  /*}}}*/
989  /*p=20, npoint=79 {{{*/
990  static IssmPDouble wgt20[]={
991  0.057256499746719, 0.001501721280705, 0.001501721280705,
992  0.001501721280705, 0.020195803723819, 0.020195803723819,
993  0.020195803723819, 0.039624016090841, 0.039624016090841,
994  0.039624016090841, 0.052739185030045, 0.052739185030045,
995  0.052739185030045, 0.053043868444611, 0.053043868444611,
996  0.053043868444611, 0.042206713977986, 0.042206713977986,
997  0.042206713977986, 0.027708365070095, 0.027708365070095,
998  0.027708365070095, 0.013333849876622, 0.013333849876622,
999  0.013333849876622,-0.001094760895106,-0.001094760895106,
1000  -0.001094760895106, 0.003033053580543, 0.003033053580543,
1001  0.003033053580543, 0.028519670065604, 0.028519670065604,
1002  0.028519670065604, 0.028519670065604, 0.028519670065604,
1003  0.028519670065604, 0.008381451951650, 0.008381451951650,
1004  0.008381451951650, 0.008381451951650, 0.008381451951650,
1005  0.008381451951650, 0.044695409202580, 0.044695409202580,
1006  0.044695409202580, 0.044695409202580, 0.044695409202580,
1007  0.044695409202580, 0.014672360101834, 0.014672360101834,
1008  0.014672360101834, 0.014672360101834, 0.014672360101834,
1009  0.014672360101834, 0.031791643800640, 0.031791643800640,
1010  0.031791643800640, 0.031791643800640, 0.031791643800640,
1011  0.031791643800640, 0.001220064691226, 0.001220064691226,
1012  0.001220064691226, 0.001220064691226, 0.001220064691226,
1013  0.001220064691226, 0.017515684095300, 0.017515684095300,
1014  0.017515684095300, 0.017515684095300, 0.017515684095300,
1015  0.017515684095300, 0.006190192638113, 0.006190192638113,
1016  0.006190192638113, 0.006190192638113, 0.006190192638113,
1017  0.006190192638113};
1018  static IssmPDouble l120[]={
1019  0.333333333333333,-0.001900928704400, 0.500950464352200,
1020  0.500950464352200, 0.023574084130543, 0.488212957934729,
1021  0.488212957934729, 0.089726636099435, 0.455136681950283,
1022  0.455136681950283, 0.196007481363421, 0.401996259318289,
1023  0.401996259318289, 0.488214180481157, 0.255892909759421,
1024  0.255892909759421, 0.647023488009788, 0.176488255995106,
1025  0.176488255995106, 0.791658289326483, 0.104170855336758,
1026  0.104170855336758, 0.893862072318140, 0.053068963840930,
1027  0.053068963840930, 0.916762569607942, 0.041618715196029,
1028  0.041618715196029, 0.976836157186356, 0.011581921406822,
1029  0.011581921406822, 0.048741583664839, 0.048741583664839,
1030  0.344855770229001, 0.606402646106160, 0.344855770229001,
1031  0.606402646106160, 0.006314115948605, 0.006314115948605,
1032  0.377843269594854, 0.615842614456541, 0.377843269594854,
1033  0.615842614456541, 0.134316520547348, 0.134316520547348,
1034  0.306635479062357, 0.559048000390295, 0.306635479062357,
1035  0.559048000390295, 0.013973893962392, 0.013973893962392,
1036  0.249419362774742, 0.736606743262866, 0.249419362774742,
1037  0.736606743262866, 0.075549132909764, 0.075549132909764,
1038  0.212775724802802, 0.711675142287434, 0.212775724802802,
1039  0.711675142287434,-0.008368153208227,-0.008368153208227,
1040  0.146965436053239, 0.861402717154987, 0.146965436053239,
1041  0.861402717154987, 0.026686063258714, 0.026686063258714,
1042  0.137726978828923, 0.835586957912363, 0.137726978828923,
1043  0.835586957912363, 0.010547719294141, 0.010547719294141,
1044  0.059696109149007, 0.929756171556853, 0.059696109149007,
1045  0.929756171556853};
1046  static IssmPDouble l220[]={
1047  0.333333333333333, 0.500950464352200,-0.001900928704400,
1048  0.500950464352200, 0.488212957934729, 0.023574084130543,
1049  0.488212957934729, 0.455136681950283, 0.089726636099435,
1050  0.455136681950283, 0.401996259318289, 0.196007481363421,
1051  0.401996259318289, 0.255892909759421, 0.488214180481157,
1052  0.255892909759421, 0.176488255995106, 0.647023488009788,
1053  0.176488255995106, 0.104170855336758, 0.791658289326483,
1054  0.104170855336758, 0.053068963840930, 0.893862072318140,
1055  0.053068963840930, 0.041618715196029, 0.916762569607942,
1056  0.041618715196029, 0.011581921406822, 0.976836157186356,
1057  0.011581921406822, 0.344855770229001, 0.606402646106160,
1058  0.048741583664839, 0.048741583664839, 0.606402646106160,
1059  0.344855770229001, 0.377843269594854, 0.615842614456541,
1060  0.006314115948605, 0.006314115948605, 0.615842614456541,
1061  0.377843269594854, 0.306635479062357, 0.559048000390295,
1062  0.134316520547348, 0.134316520547348, 0.559048000390295,
1063  0.306635479062357, 0.249419362774742, 0.736606743262866,
1064  0.013973893962392, 0.013973893962392, 0.736606743262866,
1065  0.249419362774742, 0.212775724802802, 0.711675142287434,
1066  0.075549132909764, 0.075549132909764, 0.711675142287434,
1067  0.212775724802802, 0.146965436053239, 0.861402717154987,
1068  -0.008368153208227,-0.008368153208227, 0.861402717154987,
1069  0.146965436053239, 0.137726978828923, 0.835586957912363,
1070  0.026686063258714, 0.026686063258714, 0.835586957912363,
1071  0.137726978828923, 0.059696109149007, 0.929756171556853,
1072  0.010547719294141, 0.010547719294141, 0.929756171556853,
1073  0.059696109149007};
1074  static IssmPDouble l320[]={
1075  0.333333333333333, 0.500950464352200, 0.500950464352200,
1076  -0.001900928704400, 0.488212957934729, 0.488212957934729,
1077  0.023574084130543, 0.455136681950283, 0.455136681950283,
1078  0.089726636099435, 0.401996259318289, 0.401996259318289,
1079  0.196007481363421, 0.255892909759421, 0.255892909759421,
1080  0.488214180481157, 0.176488255995106, 0.176488255995106,
1081  0.647023488009788, 0.104170855336758, 0.104170855336758,
1082  0.791658289326483, 0.053068963840930, 0.053068963840930,
1083  0.893862072318140, 0.041618715196029, 0.041618715196029,
1084  0.916762569607942, 0.011581921406822, 0.011581921406822,
1085  0.976836157186356, 0.606402646106160, 0.344855770229001,
1086  0.606402646106160, 0.344855770229001, 0.048741583664839,
1087  0.048741583664839, 0.615842614456541, 0.377843269594854,
1088  0.615842614456541, 0.377843269594854, 0.006314115948605,
1089  0.006314115948605, 0.559048000390295, 0.306635479062357,
1090  0.559048000390295, 0.306635479062357, 0.134316520547348,
1091  0.134316520547348, 0.736606743262866, 0.249419362774742,
1092  0.736606743262866, 0.249419362774742, 0.013973893962392,
1093  0.013973893962392, 0.711675142287434, 0.212775724802802,
1094  0.711675142287434, 0.212775724802802, 0.075549132909764,
1095  0.075549132909764, 0.861402717154987, 0.146965436053239,
1096  0.861402717154987, 0.146965436053239,-0.008368153208227,
1097  -0.008368153208227, 0.835586957912363, 0.137726978828923,
1098  0.835586957912363, 0.137726978828923, 0.026686063258714,
1099  0.026686063258714, 0.929756171556853, 0.059696109149007,
1100  0.929756171556853, 0.059696109149007, 0.010547719294141,
1101  0.010547719294141};
1102  /*}}}*/
1103 
1104  static IssmPDouble* wgtp[MAX_TRIA_SYM_ORD]={
1105  wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 ,
1106  wgt6 ,wgt7 ,wgt8 ,wgt9 ,wgt10,
1107  wgt11,wgt12,wgt13,wgt14,wgt15,
1108  wgt16,wgt17,wgt18,wgt19,wgt20};
1109  static IssmPDouble* l1p [MAX_TRIA_SYM_ORD]={
1110  l11 ,l12 ,l13 ,l14 ,l15 ,
1111  l16 ,l17 ,l18 ,l19 ,l110 ,
1112  l111 ,l112 ,l113 ,l114 ,l115 ,
1113  l116 ,l117 ,l118 ,l119 ,l120 };
1114  static IssmPDouble* l2p [MAX_TRIA_SYM_ORD]={
1115  l21 ,l22 ,l23 ,l24 ,l25 ,
1116  l26 ,l27 ,l28 ,l29 ,l210 ,
1117  l211 ,l212 ,l213 ,l214 ,l215 ,
1118  l216 ,l217 ,l218 ,l219 ,l220 };
1119  static IssmPDouble* l3p [MAX_TRIA_SYM_ORD]={
1120  l31 ,l32 ,l33 ,l34 ,l35 ,
1121  l36 ,l37 ,l38 ,l39 ,l310 ,
1122  l311 ,l312 ,l313 ,l314 ,l315 ,
1123  l316 ,l317 ,l318 ,l319 ,l320 };
1124 
1125  static int np[MAX_TRIA_SYM_ORD]={sizeof(wgt1 )/sizeof(IssmPDouble),
1126  sizeof(wgt2 )/sizeof(IssmPDouble),
1127  sizeof(wgt3 )/sizeof(IssmPDouble),
1128  sizeof(wgt4 )/sizeof(IssmPDouble),
1129  sizeof(wgt5 )/sizeof(IssmPDouble),
1130  sizeof(wgt6 )/sizeof(IssmPDouble),
1131  sizeof(wgt7 )/sizeof(IssmPDouble),
1132  sizeof(wgt8 )/sizeof(IssmPDouble),
1133  sizeof(wgt9 )/sizeof(IssmPDouble),
1134  sizeof(wgt10)/sizeof(IssmPDouble),
1135  sizeof(wgt11)/sizeof(IssmPDouble),
1136  sizeof(wgt12)/sizeof(IssmPDouble),
1137  sizeof(wgt13)/sizeof(IssmPDouble),
1138  sizeof(wgt14)/sizeof(IssmPDouble),
1139  sizeof(wgt15)/sizeof(IssmPDouble),
1140  sizeof(wgt16)/sizeof(IssmPDouble),
1141  sizeof(wgt17)/sizeof(IssmPDouble),
1142  sizeof(wgt18)/sizeof(IssmPDouble),
1143  sizeof(wgt19)/sizeof(IssmPDouble),
1144  sizeof(wgt20)/sizeof(IssmPDouble)};
1145 
1146  // _printf0_("GaussLegendreTria: iord=" << iord << "\n");
1147 
1148  /* check to see if Gauss points need to be calculated */
1149  if (iord <= MAX_TRIA_SYM_ORD) {
1150 
1151  /* copy the points from the static arrays (noting that the pointers
1152  could be returned directly, but then the calling function would
1153  have to know to not free them) */
1154 
1155  *pngaus=np[iord-1];
1156 
1157  *pl1 =xNew<IssmDouble>(*pngaus);
1158  *pl2 =xNew<IssmDouble>(*pngaus);
1159  *pl3 =xNew<IssmDouble>(*pngaus);
1160  *pwgt =xNew<IssmDouble>(*pngaus);
1161 
1162  for (i=0; i<*pngaus; i++) {
1163  (*pl1 )[i]=l1p [iord-1][i];
1164  (*pl2 )[i]=l2p [iord-1][i];
1165  (*pl3 )[i]=l3p [iord-1][i];
1166  (*pwgt)[i]=wgtp[iord-1][i];
1167  }
1168  }
1169  else {
1170 
1171  /* calculate the Gauss points from the collapsed quadrilateral */
1172  nigaus =iord/2+1;
1173  *pngaus=nigaus*nigaus;
1174 
1175  *pl1 =xNew<IssmDouble>(*pngaus);
1176  *pl2 =xNew<IssmDouble>(*pngaus);
1177  *pl3 =xNew<IssmDouble>(*pngaus);
1178  *pwgt =xNew<IssmDouble>(*pngaus);
1179 
1180  /* get the gauss points in each direction */
1181  GaussLegendreLinear(&xgaus, &xwgt, nigaus);
1182 
1183  egaus=xgaus;
1184  ewgt =xwgt;
1185 
1186  /* collapse the gauss points into the triangle and transform into
1187  area coordinates */
1188  ipt=0;
1189  for (j=0; j<nigaus; j++) {
1190  for (i=0; i<nigaus; i++) {
1191  xi = 1./2.*(1.-egaus[j])*xgaus[i];
1192  eta =SQRT3/2.*(1.+egaus[j]);
1193  (*pwgt)[ipt]=xwgt[i]*ewgt[j]*(SQRT3/4.*(1.-egaus[j]));
1194 
1195  (*pl1 )[ipt]=(1.-xi-eta/SQRT3)/2.;
1196  (*pl2 )[ipt]=(1.+xi-eta/SQRT3)/2.;
1197  (*pl3 )[ipt]= eta/SQRT3;
1198 
1199  ipt++;
1200  }
1201  }
1202  xDelete<IssmPDouble>(xwgt );
1203  xDelete<IssmPDouble>(xgaus);
1204  }
1205 
1206  // _printf0_("GaussLegendreTria - ngaus=" << *pngaus << "\n");
1207  // for (i=0; i<*pngaus; i++)
1208  // _printf0_("i=" << i << ": l1gaus=" << (*pl1 )[i] << ",l2gaus=" << (*pl2 )[i] << ",l3gaus=" << (*pl3 )[i]<< ",wgt=" << (*pwgt)[i]<< "\n\n");
1209 
1210  return;
1211 }/*}}}*/

◆ GaussLegendreTetra()

void GaussLegendreTetra ( int *  pngaus,
IssmDouble **  pl1,
IssmDouble **  pl2,
IssmDouble **  pl3,
IssmDouble **  pl4,
IssmDouble **  pwgt,
int  iord 
)

Definition at line 1212 of file GaussPoints.cpp.

1212  {/*{{{*/
1213  /* Gauss quadrature points for the tetrahedron.
1214 
1215  p=2-3 points from Y. Jinyun, "Symmetric Gaussian Quadrature
1216  Formulae for Tetrahedronal Regions", Computer Methods in Applied
1217  Mechanics and Engineering, Vol. 43, pp. 349-353 (1984).
1218 
1219  p=4-6 points from P. Keast, "Moderate-Degree Tetrahedral
1220  Quadrature Formulas", Computer Methods in Applied Mechanics and
1221  Engineering, Vol. 55, pp. 339-348 (1986).*/
1222  _assert_(iord>0);
1223 
1224  /*Intermediaries*/
1225  int i,j,k,ipt,nigaus;
1226  IssmPDouble xi,eta,zeta;
1227  IssmPDouble *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt;
1228 
1229  /*Hardcoded Gauss points definition*/
1230  /*p= 1, npoint= 1 {{{*/
1231  static IssmPDouble wgt1[]={
1232  1.000000000000000};
1233  static IssmPDouble l11[]={
1234  0.250000000000000};
1235  static IssmPDouble l21[]={
1236  0.250000000000000};
1237  static IssmPDouble l31[]={
1238  0.250000000000000};
1239  static IssmPDouble l41[]={
1240  0.250000000000000};
1241  /*}}}*/
1242  /*p= 2, npoint= 4 {{{*/
1243 
1244  static IssmPDouble wgt2[]={
1245  0.250000000000000, 0.250000000000000, 0.250000000000000,
1246  0.250000000000000};
1247  static IssmPDouble l12[]={
1248  0.585410196624969, 0.138196601125011, 0.138196601125011,
1249  0.138196601125011};
1250  static IssmPDouble l22[]={
1251  0.138196601125011, 0.585410196624969, 0.138196601125011,
1252  0.138196601125011};
1253  static IssmPDouble l32[]={
1254  0.138196601125011, 0.138196601125011, 0.585410196624969,
1255  0.138196601125011};
1256  static IssmPDouble l42[]={
1257  0.138196601125011, 0.138196601125011, 0.138196601125011,
1258  0.585410196624969};
1259  /*}}}*/
1260  /*p= 3, npoint= 5 {{{*/
1261  static IssmPDouble wgt3[]={
1262  -0.800000000000000, 0.450000000000000, 0.450000000000000,
1263  0.450000000000000, 0.450000000000000};
1264  static IssmPDouble l13[]={
1265  0.250000000000000, 0.500000000000000, 0.166666666666667,
1266  0.166666666666667, 0.166666666666667};
1267  static IssmPDouble l23[]={
1268  0.250000000000000, 0.166666666666667, 0.500000000000000,
1269  0.166666666666667, 0.166666666666667};
1270  static IssmPDouble l33[]={
1271  0.250000000000000, 0.166666666666667, 0.166666666666667,
1272  0.500000000000000, 0.166666666666667};
1273  static IssmPDouble l43[]={
1274  0.250000000000000, 0.166666666666667, 0.166666666666667,
1275  0.166666666666667, 0.500000000000000};
1276  /*}}}*/
1277  /*p= 4, npoint=11 {{{*/
1278 
1279  static IssmPDouble wgt4[]={
1280  -0.013155555555556, 0.007622222222222, 0.007622222222222,
1281  0.007622222222222, 0.007622222222222, 0.024888888888889,
1282  0.024888888888889, 0.024888888888889, 0.024888888888889,
1283  0.024888888888889, 0.024888888888889};
1284  static IssmPDouble l14[]={
1285  0.250000000000000, 0.785714285714286, 0.071428571428571,
1286  0.071428571428571, 0.071428571428571, 0.399403576166799,
1287  0.399403576166799, 0.399403576166799, 0.100596423833201,
1288  0.100596423833201, 0.100596423833201};
1289  static IssmPDouble l24[]={
1290  0.250000000000000, 0.071428571428571, 0.785714285714286,
1291  0.071428571428571, 0.071428571428571, 0.399403576166799,
1292  0.100596423833201, 0.100596423833201, 0.399403576166799,
1293  0.399403576166799, 0.100596423833201};
1294  static IssmPDouble l34[]={
1295  0.250000000000000, 0.071428571428571, 0.071428571428571,
1296  0.785714285714286, 0.071428571428571, 0.100596423833201,
1297  0.399403576166799, 0.100596423833201, 0.399403576166799,
1298  0.100596423833201, 0.399403576166799};
1299  static IssmPDouble l44[]={
1300  0.250000000000000, 0.071428571428571, 0.071428571428571,
1301  0.071428571428571, 0.785714285714286, 0.100596423833201,
1302  0.100596423833201, 0.399403576166799, 0.100596423833201,
1303  0.399403576166799, 0.399403576166799};
1304  /*}}}*/
1305  /*p= 5, npoint=15 {{{*/
1306 
1307  static IssmPDouble wgt5[]={
1308  0.030283678097089, 0.006026785714286, 0.006026785714286,
1309  0.006026785714286, 0.006026785714286, 0.011645249086029,
1310  0.011645249086029, 0.011645249086029, 0.011645249086029,
1311  0.010949141561386, 0.010949141561386, 0.010949141561386,
1312  0.010949141561386, 0.010949141561386, 0.010949141561386};
1313  static IssmPDouble l15[]={
1314  0.250000000000000, 0.000000000000000, 0.333333333333333,
1315  0.333333333333333, 0.333333333333333, 0.727272727272727,
1316  0.090909090909091, 0.090909090909091, 0.090909090909091,
1317  0.066550153573664, 0.066550153573664, 0.066550153573664,
1318  0.433449846426336, 0.433449846426336, 0.433449846426336};
1319  static IssmPDouble l25[]={
1320  0.250000000000000, 0.333333333333333, 0.000000000000000,
1321  0.333333333333333, 0.333333333333333, 0.090909090909091,
1322  0.727272727272727, 0.090909090909091, 0.090909090909091,
1323  0.066550153573664, 0.433449846426336, 0.433449846426336,
1324  0.066550153573664, 0.066550153573664, 0.433449846426336};
1325  static IssmPDouble l35[]={
1326  0.250000000000000, 0.333333333333333, 0.333333333333333,
1327  0.000000000000000, 0.333333333333333, 0.090909090909091,
1328  0.090909090909091, 0.727272727272727, 0.090909090909091,
1329  0.433449846426336, 0.066550153573664, 0.433449846426336,
1330  0.066550153573664, 0.433449846426336, 0.066550153573664};
1331  static IssmPDouble l45[]={
1332  0.250000000000000, 0.333333333333333, 0.333333333333333,
1333  0.333333333333333, 0.000000000000000, 0.090909090909091,
1334  0.090909090909091, 0.090909090909091, 0.727272727272727,
1335  0.433449846426336, 0.433449846426336, 0.066550153573664,
1336  0.433449846426336, 0.066550153573664, 0.066550153573664};
1337  /*}}}*/
1338  /*p= 6, npoint=24 {{{*/
1339 
1340  static IssmPDouble wgt6[]={
1341  0.006653791709695, 0.006653791709695, 0.006653791709695,
1342  0.006653791709695, 0.001679535175887, 0.001679535175887,
1343  0.001679535175887, 0.001679535175887, 0.009226196923942,
1344  0.009226196923942, 0.009226196923942, 0.009226196923942,
1345  0.008035714285714, 0.008035714285714, 0.008035714285714,
1346  0.008035714285714, 0.008035714285714, 0.008035714285714,
1347  0.008035714285714, 0.008035714285714, 0.008035714285714,
1348  0.008035714285714, 0.008035714285714, 0.008035714285714};
1349  static IssmPDouble l16[]={
1350  0.356191386222545, 0.214602871259152, 0.214602871259152,
1351  0.214602871259152, 0.877978124396166, 0.040673958534611,
1352  0.040673958534611, 0.040673958534611, 0.032986329573173,
1353  0.322337890142276, 0.322337890142276, 0.322337890142276,
1354 
1355  0.063661001875018, 0.063661001875018, 0.063661001875018,
1356  0.063661001875018, 0.063661001875018, 0.063661001875018,
1357  0.269672331458316, 0.603005664791649, 0.269672331458316,
1358  0.603005664791649, 0.269672331458316, 0.603005664791649};
1359  static IssmPDouble l26[]={
1360  0.214602871259152, 0.356191386222545, 0.214602871259152,
1361  0.214602871259152, 0.040673958534611, 0.877978124396166,
1362  0.040673958534611, 0.040673958534611, 0.322337890142276,
1363  0.032986329573173, 0.322337890142276, 0.322337890142276,
1364 
1365  0.063661001875018, 0.063661001875018, 0.269672331458316,
1366  0.603005664791649, 0.269672331458316, 0.603005664791649,
1367  0.063661001875018, 0.063661001875018, 0.063661001875018,
1368  0.063661001875018, 0.603005664791649, 0.269672331458316};
1369  static IssmPDouble l36[]={
1370  0.214602871259152, 0.214602871259152, 0.356191386222545,
1371  0.214602871259152, 0.040673958534611, 0.040673958534611,
1372  0.877978124396166, 0.040673958534611, 0.322337890142276,
1373  0.322337890142276, 0.032986329573173, 0.322337890142276,
1374 
1375  0.269672331458316, 0.603005664791649, 0.063661001875018,
1376  0.063661001875018, 0.603005664791649, 0.269672331458316,
1377  0.063661001875018, 0.063661001875018, 0.603005664791649,
1378  0.269672331458316, 0.063661001875018, 0.063661001875018};
1379  static IssmPDouble l46[]={
1380  0.214602871259152, 0.214602871259152, 0.214602871259152,
1381  0.356191386222545, 0.040673958534611, 0.040673958534611,
1382  0.040673958534611, 0.877978124396166, 0.322337890142276,
1383  0.322337890142276, 0.322337890142276, 0.032986329573173,
1384 
1385  0.603005664791649, 0.269672331458316, 0.603005664791649,
1386  0.269672331458316, 0.063661001875018, 0.063661001875018,
1387  0.603005664791649, 0.269672331458316, 0.063661001875018,
1388  0.063661001875018, 0.063661001875018, 0.063661001875018};
1389  /*}}}*/
1390 
1391  static IssmPDouble* wgtp[MAX_TETRA_SYM_ORD]={wgt1,wgt2,wgt3,wgt4,wgt5,wgt6};
1392  static IssmPDouble* l1p [MAX_TETRA_SYM_ORD]={l11 ,l12 ,l13 ,l14 ,l15 ,l16 };
1393  static IssmPDouble* l2p [MAX_TETRA_SYM_ORD]={l21 ,l22 ,l23 ,l24 ,l25 ,l26 };
1394  static IssmPDouble* l3p [MAX_TETRA_SYM_ORD]={l31 ,l32 ,l33 ,l34 ,l35 ,l36 };
1395  static IssmPDouble* l4p [MAX_TETRA_SYM_ORD]={l41 ,l42 ,l43 ,l44 ,l45 ,l46 };
1396 
1397  static int np[MAX_TETRA_SYM_ORD]={sizeof(wgt1 )/sizeof(IssmPDouble),
1398  sizeof(wgt2 )/sizeof(IssmPDouble),
1399  sizeof(wgt3 )/sizeof(IssmPDouble),
1400  sizeof(wgt4 )/sizeof(IssmPDouble),
1401  sizeof(wgt5 )/sizeof(IssmPDouble),
1402  sizeof(wgt6 )/sizeof(IssmPDouble)};
1403 
1404  // _printf0_("GaussLegendreTetra: iord=" << iord << "\n");
1405 
1406  /* check to see if Gauss points need to be calculated */
1407  if (iord <= MAX_TETRA_SYM_ORD) {
1408 
1409  /* copy the points from the static arrays (noting that the pointers
1410  could be returned directly, but then the calling function would
1411  have to know to not free them), and multiply the weights by the
1412  volume of the parametric tetrahedron */
1413 
1414  *pngaus=np[iord-1];
1415 
1416  *pl1 =xNew<IssmDouble>(*pngaus);
1417  *pl2 =xNew<IssmDouble>(*pngaus);
1418  *pl3 =xNew<IssmDouble>(*pngaus);
1419  *pl4 =xNew<IssmDouble>(*pngaus);
1420  *pwgt =xNew<IssmDouble>(*pngaus);
1421 
1422  for (i=0; i<*pngaus; i++) {
1423  (*pl1 )[i]=l1p [iord-1][i];
1424  (*pl2 )[i]=l2p [iord-1][i];
1425  (*pl3 )[i]=l3p [iord-1][i];
1426  (*pl4 )[i]=l4p [iord-1][i];
1427  (*pwgt)[i]=wgtp[iord-1][i];
1428  }
1429  }
1430  else {
1431 
1432  /* calculate the Gauss points from the collapsed hexahedron */
1433  nigaus =iord/2+1;
1434  *pngaus=nigaus*nigaus*nigaus;
1435 
1436  *pl1 =xNew<IssmDouble>(*pngaus);
1437  *pl2 =xNew<IssmDouble>(*pngaus);
1438  *pl3 =xNew<IssmDouble>(*pngaus);
1439  *pl4 =xNew<IssmDouble>(*pngaus);
1440  *pwgt =xNew<IssmDouble>(*pngaus);
1441 
1442  /* get the gauss points in each direction */
1443  GaussLegendreLinear(&xgaus, &xwgt, nigaus);
1444 
1445  egaus=xgaus;
1446  ewgt =xwgt;
1447  zgaus=xgaus;
1448  zwgt =xwgt;
1449 
1450  /* collapse the gauss points into the tetrahedron and transform into
1451  volume coordinates */
1452  ipt=0;
1453  for (k=0; k<nigaus; k++) {
1454  for (j=0; j<nigaus; j++) {
1455  for (i=0; i<nigaus; i++) {
1456  xi =1./4.*(1.-egaus[j])*(1.-zgaus[k])*xgaus[i];
1457  eta =1./4./SQRT3
1458  *(5.+3.*egaus[j]-zgaus[k]-3.*egaus[j]*zgaus[k]);
1459  zeta =sqrt(2./3.)*(1.+zgaus[k]);
1460  (*pwgt)[ipt]=xwgt[i]*ewgt[j]*zwgt[k]*(SQRT2/16.*(1.-egaus[j])*pow(1.-zgaus[k],2));
1461 
1462  (*pl1 )[ipt]=(1.-xi-eta/SQRT3-zeta/sqrt(6.))/2.;
1463  (*pl2 )[ipt]=(1.+xi-eta/SQRT3-zeta/sqrt(6.))/2.;
1464  (*pl3 )[ipt]=( eta -zeta/sqrt(8.))/SQRT3;
1465  (*pl4 )[ipt]=( zeta/sqrt(8.))*SQRT3;
1466 
1467  ipt++;
1468  }
1469  }
1470  }
1471  xDelete<IssmPDouble>(xwgt );
1472  xDelete<IssmPDouble>(xgaus);
1473  }
1474 }/*}}}*/

◆ GaussLobatto()

void GaussLobatto ( IssmPDouble **  pxgaus,
IssmPDouble **  pxwgt,
int  ngaus 
)

Definition at line 1475 of file GaussPoints.cpp.

1475  {/*{{{*/
1476  /*Gauss-Lobatto quadrature points.
1477 
1478  The recurrence coefficients for Legendre polynomials on (-1,1)
1479  are defined (from the ORTHPOL subroutine RECUR with ipoly=1) as:
1480 
1481  alpha(i)=0.
1482  beta (i)=1./(4.-1./(i-1)^2))
1483 
1484  and then modified for the Gauss-Lobatto quadrature rule on (-1,1)
1485  (from the ORTHPOL subroutine LOB).
1486 
1487  For degree p, the required number of Gauss-Lobatto points is
1488  n>=(p+1)/2+1 (one more than Gauss-Legendre).*/
1489  _assert_(ngaus>0);
1490 
1491  int i;
1492  IssmPDouble *alpha,*beta;
1493  IssmPDouble left=-1.,right= 1.;
1494  IssmPDouble p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;
1495 
1496  /*p= 1, npoint= 1 (Gauss-Legendre)*/
1497  static IssmPDouble wgt1[]={2.000000000000000};
1498  static IssmPDouble xi1[]={0.000000000000000};
1499 
1500  /*p= 1, npoint= 2*/
1501  static IssmPDouble wgt2[]={1.000000000000000, 1.000000000000000};
1502  static IssmPDouble xi2[]={-1.000000000000000, 1.000000000000000};
1503 
1504  /*p= 3, npoint= 3*/
1505  static IssmPDouble wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
1506  static IssmPDouble xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
1507 
1508  /*p= 5, npoint= 4*/
1509  static IssmPDouble wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};
1510  static IssmPDouble xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};
1511 
1512  /*p= 7, npoint= 5*/
1513  static IssmPDouble wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};
1514  static IssmPDouble xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};
1515 
1516  static IssmPDouble* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };
1517  static IssmPDouble* xip [MAX_LINE_GLOB_PTS]={xi1 ,xi2 ,xi3 ,xi4 ,xi5 };
1518 
1519  *pxgaus =xNew<IssmPDouble>(ngaus);
1520  *pxwgt =xNew<IssmPDouble>(ngaus);
1521 
1522  /* check to see if Gauss points need to be calculated */
1523  if (ngaus <= MAX_LINE_GLOB_PTS) {
1524 
1525  /* copy the points from the static arrays (noting that the pointers
1526  could be returned directly, but then the calling function would
1527  have to know to not free them) */
1528  for (i=0; i<ngaus; i++) {
1529  (*pxgaus)[i]=xip [ngaus-1][i];
1530  (*pxwgt )[i]=wgtp[ngaus-1][i];
1531  }
1532  }
1533  else {
1534 
1535  /* calculate the Gauss points using recurrence relations */
1536  alpha=xNew<IssmPDouble>(ngaus);
1537  beta =xNew<IssmPDouble>(ngaus);
1538 
1539  /* calculate the Legendre recurrence coefficients */
1540  alpha[0]=0.;
1541  beta [0]=2.;
1542 
1543  for (i=1; i<ngaus; i++) {
1544  alpha[i]=0.;
1545  beta [i]=1./(4.-1./(i*i));
1546  }
1547 
1548  /* calculate the Gauss-Lobatto quadrature formula */
1549  for (i=0; i<ngaus-1; i++) {
1550  pm1l=p0l;
1551  p0l=p1l;
1552  pm1r=p0r;
1553  p0r=p1r;
1554  p1l=(left -alpha[i])*p0l-beta[i]*pm1l;
1555  p1r=(right-alpha[i])*p0r-beta[i]*pm1r;
1556  }
1557 
1558  /* Normalize system to prevent underflow:
1559  [ p1l p0l ]{ a } = {left *p1l}
1560  [ p1r p0r ]{ b } {right*p1r}
1561  dividing by p1l in the first equation and p1r in the second. */
1562 
1563  // det=p1l*p0r-p1r*p0l;
1564  det=p0r/p1r-p0l/p1l;
1565  // alpha[ngaus-1]=(left*p1l*p0r-right*p1r*p0l)/det;
1566  // beta [ngaus-1]=(right-left)*p1l*p1r/det;
1567  alpha[ngaus-1]=(left *(p0r/p1r)-right*(p0l/p1l))/det;
1568  beta [ngaus-1]=(right -left )/det;
1569 
1570  /* calculate the Gauss points */
1571  GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
1572  xDelete<IssmPDouble>(beta );
1573  xDelete<IssmPDouble>(alpha);
1574  }
1575 
1576 }/*}}}*/

◆ GaussRecur()

void GaussRecur ( IssmPDouble *  zero,
IssmPDouble *  weight,
int  n,
IssmPDouble *  alpha,
IssmPDouble *  beta 
)

Definition at line 1577 of file GaussPoints.cpp.

1577  {/*{{{*/
1578  /*Gauss quadrature points from recursion coefficients.
1579  *
1580  *The routine uses the algorithm from the ORTHPOL routine GAUSS, which
1581  *finds the eigenvalues of a tridiagonal matrix.*/
1582 
1583  /*Intermediaries*/
1584  int i,j,k,l,m,ii,mml,iter;
1585  IssmPDouble p,g,r,s,c,f,b;
1586  IssmPDouble* work;
1587 
1588  if (n==1){
1589  zero[0] =alpha[0];
1590  weight[0]=beta[0];
1591  return;
1592  }
1593 
1594  work=xNew<IssmPDouble>(n);
1595 
1596  zero[0] =alpha[0];
1597  weight[0]=1.;
1598  work[n-1]=0.;
1599  for (i=1; i<=n-1; i++){
1600  zero[i]=alpha[i];
1601  work[i-1]=sqrt(beta[i]);
1602  weight[i]=0;
1603  }
1604 
1605  for (l=0; l<=n-1; l++){
1606  iter=0;
1607  do {
1608 
1609  /* Look for a small subdiagonal element. */
1610  for (m=l; m<=n-1; m++) {
1611  if (m == n-1) break;
1612  if (fabs(work[m])
1613  <= DBL_EPSILON*(fabs(zero[m])+fabs(zero[m+1])))
1614  break;
1615  }
1616  p=zero[l];
1617  if (m==l) break;
1618  ++iter;
1619 
1620  /* Form shift. */
1621  g=(zero[l+1]-p)/(2.*work[l]);
1622  r=sqrt(g*g+1.);
1623  // g=zero[m]-p+work[l]/(g+FortranSign(r,g));
1624  g=zero[m]-p+work[l]/(g+(g>=0 ? fabs(r) : -fabs(r)));
1625  s=1.;
1626  c=1.;
1627  p=0.;
1628  mml=m-l;
1629 
1630  /* For i=m-1 step -1 until l do ... */
1631  for (ii=1; ii<=mml; ii++) {
1632  i=m-ii;
1633  f=s*work[i];
1634  b=c*work[i];
1635  if (fabs(f) >= fabs(g)) {
1636  c=g/f;
1637  r=sqrt(c*c+1.);
1638  work[i+1]=f*r;
1639  s=1./r;
1640  c*=s;
1641  }
1642  else {
1643  s=f/g;
1644  r=sqrt(s*s+1.);
1645  work[i+1]=g*r;
1646  c=1./r;
1647  s*=c;
1648  }
1649  g=zero[i+1]-p;
1650  r=(zero[i]-g)*s+2.*c*b;
1651  p=s*r;
1652  zero[i+1]=g+p;
1653  g=c*r-b;
1654 
1655  /* Form first component of vector. */
1656  f=weight[i+1];
1657  weight[i+1]=s*weight[i]+c*f;
1658  weight[i ]=c*weight[i]-s*f;
1659  }
1660  zero[l]-=p;
1661  work[l]=g;
1662  work[m]=0.;
1663  } while (iter < MAX_GAUS_ITER);
1664  if (iter >= MAX_GAUS_ITER) {
1665  xDelete<IssmPDouble>(work);
1666  _error_("Max iterations exceeded for l=" << MAX_GAUS_ITER);
1667  }
1668  }
1669 
1670  /* Order eigenvalues and eigenvectors. */
1671  for (i=0;i<n-1;i++) {
1672  k=i;
1673  p=zero[i];
1674  for (j=i+1;j<n;j++){
1675  if (zero[j] < p){
1676  k=j;
1677  p=zero[j];
1678  }
1679  }
1680  if (k > i){
1681  p=zero[i];
1682  zero[i]=zero[k];
1683  zero[k]=p;
1684  p=weight[i];
1685  weight[i]=weight[k];
1686  weight[k]=p;
1687  }
1688  }
1689  for (i=0;i<n;i++){
1690  weight[i]=beta[0]*weight[i]*weight[i];
1691  }
1692 
1693  /*Cleanup*/
1694  xDelete<IssmPDouble>(work);
1695 
1696 }/*}}}*/
GaussRecur
void GaussRecur(IssmPDouble *zero, IssmPDouble *weight, int n, IssmPDouble *alpha, IssmPDouble *beta)
Definition: GaussPoints.cpp:1577
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
MAX_TRIA_SYM_ORD
#define MAX_TRIA_SYM_ORD
Definition: GaussPoints.h:11
MAX_LINE_GLOB_PTS
#define MAX_LINE_GLOB_PTS
Definition: GaussPoints.h:15
MAX_LINE_GAUS_PTS
#define MAX_LINE_GAUS_PTS
Definition: GaussPoints.h:9
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
bamg::det
long long det(const I2 &a, const I2 &b, const I2 &c)
Definition: det.h:8
GaussLegendreLinear
void GaussLegendreLinear(IssmPDouble **pxgaus, IssmPDouble **pxwgt, int ngaus)
Definition: GaussPoints.cpp:11
MAX_GAUS_ITER
#define MAX_GAUS_ITER
Definition: GaussPoints.h:17
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
SQRT2
#define SQRT2
Definition: constants.h:9
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
SQRT3
#define SQRT3
Definition: constants.h:10
MAX_TETRA_SYM_ORD
#define MAX_TETRA_SYM_ORD
Definition: GaussPoints.h:13
Generated on Thu Jul 2 2020 08:09:25 for Ice Sheet System Model by   doxygen 1.8.19