Ice Sheet System Model  4.18
Code documentation
Public Member Functions | Data Fields | Private Attributes
GaussPenta Class Reference

#include <GaussPenta.h>

Inheritance diagram for GaussPenta:
Gauss

Public Member Functions

 GaussPenta ()
 
 GaussPenta (int order_horiz, int order_vert)
 
 GaussPenta (int index1, int index2, int order)
 
 GaussPenta (int index1, int index2, int index3, int order)
 
 GaussPenta (int index1, int index2, int index3, int index4, int order_horiz, int order_vert)
 
 GaussPenta (int index, IssmDouble r1, IssmDouble r2, bool maintlyfloating, int order)
 
 GaussPenta (IssmDouble area_coordinates[4][3], int order_horiz, int order_vert)
 
 GaussPenta (IssmDouble area_coordinates[2][3], int order_horiz)
 
 ~GaussPenta ()
 
int begin (void)
 
void Echo (void)
 
int end (void)
 
int Enum (void)
 
void GaussFaceTria (int index1, int index2, int index3, int order)
 
void GaussNode (int finitelement, int iv)
 
void GaussPoint (int ig)
 
void GaussVertex (int iv)
 
void SynchronizeGaussBase (Gauss *gauss)
 
- Public Member Functions inherited from Gauss
virtual ~Gauss ()
 

Data Fields

IssmDouble coord1
 
IssmDouble coord2
 
IssmDouble coord3
 
IssmDouble coord4
 
- Data Fields inherited from Gauss
IssmDouble weight
 

Private Attributes

int numgauss
 
IssmDouble * weights
 
IssmDouble * coords1
 
IssmDouble * coords2
 
IssmDouble * coords3
 
IssmDouble * coords4
 

Detailed Description

Definition at line 13 of file GaussPenta.h.

Constructor & Destructor Documentation

◆ GaussPenta() [1/8]

GaussPenta::GaussPenta ( )

Definition at line 16 of file GaussPenta.cpp.

16  {/*{{{*/
17 
18  numgauss=-1;
19 
20  weights=NULL;
21  coords1=NULL;
22  coords2=NULL;
23  coords3=NULL;
24  coords4=NULL;
25 
26  weight=UNDEF;
27  coord1=UNDEF;
28  coord2=UNDEF;
29  coord3=UNDEF;
30  coord4=UNDEF;
31 }

◆ GaussPenta() [2/8]

GaussPenta::GaussPenta ( int  order_horiz,
int  order_vert 
)

Definition at line 33 of file GaussPenta.cpp.

33  {/*{{{*/
34 
35  /*Intermediaries*/
36  int ighoriz,igvert;
37  int numgauss_horiz;
38  int numgauss_vert;
39  IssmDouble *coords1_horiz = NULL;
40  IssmDouble *coords2_horiz = NULL;
41  IssmDouble *coords3_horiz = NULL;
42  IssmDouble *weights_horiz = NULL;
43  double *coords_vert = NULL;
44  double *weights_vert = NULL;
45 
46  /*Get gauss points*/
47  GaussLegendreTria(&numgauss_horiz,&coords1_horiz,&coords2_horiz,&coords3_horiz,&weights_horiz,order_horiz);
48  GaussLegendreLinear(&coords_vert,&weights_vert,order_vert);
49  numgauss_vert=order_vert;
50 
51  /*Allocate GaussPenta fields*/
52  numgauss=numgauss_horiz*numgauss_vert;
53  coords1=xNew<IssmDouble>(numgauss);
54  coords2=xNew<IssmDouble>(numgauss);
55  coords3=xNew<IssmDouble>(numgauss);
56  coords4=xNew<IssmDouble>(numgauss);
57  weights=xNew<IssmDouble>(numgauss);
58 
59  /*Combine Horizontal and vertical points*/
60  for (ighoriz=0; ighoriz<numgauss_horiz; ighoriz++){
61  for (igvert=0; igvert<numgauss_vert; igvert++){
62  coords1[numgauss_vert*ighoriz+igvert]=coords1_horiz[ighoriz];
63  coords2[numgauss_vert*ighoriz+igvert]=coords2_horiz[ighoriz];
64  coords3[numgauss_vert*ighoriz+igvert]=coords3_horiz[ighoriz];
65  coords4[numgauss_vert*ighoriz+igvert]=coords_vert[igvert];
66  weights[numgauss_vert*ighoriz+igvert]=weights_horiz[ighoriz]*weights_vert[igvert];
67  }
68  }
69 
70  /*Initialize static fields as undefinite*/
71  weight=UNDEF;
72  coord1=UNDEF;
73  coord2=UNDEF;
74  coord3=UNDEF;
75  coord4=UNDEF;
76 
77  /*Clean up*/
78  xDelete<IssmDouble>(coords1_horiz);
79  xDelete<IssmDouble>(coords2_horiz);
80  xDelete<IssmDouble>(coords3_horiz);
81  xDelete<double>(coords_vert);
82  xDelete<IssmDouble>(weights_horiz);
83  xDelete<double>(weights_vert);
84 }

◆ GaussPenta() [3/8]

GaussPenta::GaussPenta ( int  index1,
int  index2,
int  order 
)

Definition at line 86 of file GaussPenta.cpp.

86  {/*{{{*/
87 
88  /*Intermediaties*/
89  double *seg_coords = NULL;
90  double *seg_weights = NULL;
91  int i;
92 
93  /*Get Segment gauss points*/
94  numgauss=order;
95  GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);
96 
97  /*Allocate GaussPenta fields*/
98  coords1=xNew<IssmDouble>(numgauss);
99  coords2=xNew<IssmDouble>(numgauss);
100  coords3=xNew<IssmDouble>(numgauss);
101  coords4=xNew<IssmDouble>(numgauss);
102  weights=xNew<IssmDouble>(numgauss);
103 
104  if(index1==0 && index2==3){
105  for(i=0;i<numgauss;i++) coords1[i]=1.0;
106  for(i=0;i<numgauss;i++) coords2[i]=0.0;
107  for(i=0;i<numgauss;i++) coords3[i]=0.0;
108  for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
109  for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
110  }
111  else if (index1==1 && index2==4){
112  for(i=0;i<numgauss;i++) coords1[i]=0.0;
113  for(i=0;i<numgauss;i++) coords2[i]=1.0;
114  for(i=0;i<numgauss;i++) coords3[i]=0.0;
115  for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
116  for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
117  }
118  else if (index1==2 && index2==5){
119  for(i=0;i<numgauss;i++) coords1[i]=0.0;
120  for(i=0;i<numgauss;i++) coords2[i]=0.0;
121  for(i=0;i<numgauss;i++) coords3[i]=1.0;
122  for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
123  for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
124  }
125  else{
126  _error_("Penta not supported yet");
127  }
128 
129  /*Initialize static fields as undefined*/
130  weight=UNDEF;
131  coord1=UNDEF;
132  coord2=UNDEF;
133  coord3=UNDEF;
134  coord4=UNDEF;
135 
136  /*clean up*/
137  xDelete<double>(seg_coords);
138  xDelete<double>(seg_weights);
139 
140 }

◆ GaussPenta() [4/8]

GaussPenta::GaussPenta ( int  index1,
int  index2,
int  index3,
int  order 
)

Definition at line 142 of file GaussPenta.cpp.

142  {/*{{{*/
143 
144  /*Basal Tria*/
145  if(index1==0 && index2==1 && index3==2){
146 
147  /*Get GaussTria*/
148  GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
149 
150  /*compute z coordinate*/
151  coords4=xNew<IssmDouble>(numgauss);
152  for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
153  }
154  /*Upper surface Tria*/
155  else if(index1==3 && index2==4 && index3==5){
156 
157  /*Get GaussTria*/
158  GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
159 
160  /*compute z coordinate*/
161  coords4=xNew<IssmDouble>(numgauss);
162  for(int i=0;i<numgauss;i++) coords4[i]=1.0;
163  }
164  else{
165  _error_("Tria not supported yet");
166  }
167 
168 }

◆ GaussPenta() [5/8]

GaussPenta::GaussPenta ( int  index1,
int  index2,
int  index3,
int  index4,
int  order_horiz,
int  order_vert 
)

Definition at line 170 of file GaussPenta.cpp.

170  {/*{{{*/
171 
172  /*Intermediaties*/
173  double *seg_horiz_coords = NULL;
174  double *seg_horiz_weights = NULL;
175  double *seg_vert_coords = NULL;
176  double *seg_vert_weights = NULL;
177  int i,j;
178 
179  /*get the gauss points using the product of two line rules*/
180  GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
181  GaussLegendreLinear(&seg_vert_coords, &seg_vert_weights, order_vert);
182 
183  /*Allocate GaussPenta fields*/
184  numgauss=order_horiz*order_vert;
185  coords1=xNew<IssmDouble>(numgauss);
186  coords2=xNew<IssmDouble>(numgauss);
187  coords3=xNew<IssmDouble>(numgauss);
188  coords4=xNew<IssmDouble>(numgauss);
189  weights=xNew<IssmDouble>(numgauss);
190 
191  /*Quads: get the gauss points using the product of two line rules */
192  if(index1==0 && index2==1 && index3==4 && index4==3){
193  for(i=0;i<order_horiz;i++){
194  for(j=0;j<order_vert;j++){
195  coords1[i*order_vert+j]= 0.5*(1-seg_horiz_coords[i]);
196  coords2[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
197  coords3[i*order_vert+j]=0.0;
198  coords4[i*order_vert+j]=seg_vert_coords[j];
199  weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
200  }
201  }
202  }
203  else if(index1==1 && index2==2 && index3==5 && index4==4){
204  for(i=0;i<order_horiz;i++){
205  for(j=0;j<order_vert;j++){
206  coords1[i*order_vert+j]=0.0;
207  coords2[i*order_vert+j]= 0.5*(1-seg_horiz_coords[i]);
208  coords3[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
209  coords4[i*order_vert+j]=seg_vert_coords[j];
210  weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
211  }
212  }
213  }
214  else if(index1==2 && index2==0 && index3==3 && index4==5){
215  for(i=0;i<order_horiz;i++){
216  for(j=0;j<order_vert;j++){
217  coords1[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
218  coords2[i*order_vert+j]=0.0;
219  coords3[i*order_vert+j]= 0.5*(1-seg_horiz_coords[i]);
220  coords4[i*order_vert+j]=seg_vert_coords[j];
221  weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
222  }
223  }
224  }
225  else{
226  _error_("Tria not supported yet (user provided indices " << index1 << " " << index2 << " " << index3 << " " << index4 << ")");
227  }
228 
229  /*clean-up*/
230  xDelete<double>(seg_horiz_coords);
231  xDelete<double>(seg_horiz_weights);
232  xDelete<double>(seg_vert_coords);
233  xDelete<double>(seg_vert_weights);
234 }

◆ GaussPenta() [6/8]

GaussPenta::GaussPenta ( int  index,
IssmDouble  r1,
IssmDouble  r2,
bool  maintlyfloating,
int  order 
)

Definition at line 236 of file GaussPenta.cpp.

236  {/*{{{*/
237 
238  /*
239  * ^
240  * |
241  * 1|\
242  * | \
243  * | \
244  * | \
245  * | \
246  * | \
247  * | +(x,y) \
248  * | \
249  * +---------------+-->
250  * 0 1
251  *
252  */
253  int ig;
254  IssmDouble x,y;
255  IssmDouble xy_list[3][2];
256 
257  if(mainlyfloating){
258  /*Get gauss points*/
259  GaussLegendreTria(&this->numgauss,&this->coords1,&this->coords2,&this->coords3,&this->weights,order);
260 
261  xy_list[0][0]=0; xy_list[0][1]=0;
262  xy_list[1][0]=r1; xy_list[1][1]=0;
263  xy_list[2][0]=0; xy_list[2][1]=r2;
264 
265  for(ig=0;ig<this->numgauss;ig++){
266  x = this->coords1[ig]*xy_list[0][0] + this->coords2[ig]*xy_list[1][0] + this->coords3[ig]*xy_list[2][0];
267  y = this->coords1[ig]*xy_list[0][1] + this->coords2[ig]*xy_list[1][1] + this->coords3[ig]*xy_list[2][1];
268 
269  switch(index){
270  case 0:
271  this->coords1[ig] = 1.-x-y;
272  this->coords2[ig] = x;
273  this->coords3[ig] = y;
274  break;
275  case 1:
276  this->coords1[ig] = y;
277  this->coords2[ig] = 1.-x-y;
278  this->coords3[ig] = x;
279  break;
280  case 2:
281  this->coords1[ig] = x;
282  this->coords2[ig] = y;
283  this->coords3[ig] = 1.-x-y;
284  break;
285  default:
286  _error_("index "<<index<<" not supported yet");
287  }
288  this->weights[ig] = this->weights[ig]*r1*r2;
289  }
290  this->coords4=xNew<IssmDouble>(numgauss);
291  for(ig=0;ig<numgauss;ig++) this->coords4[ig]=-1.0;
292  }
293  else{
294  /*Double number of gauss points*/
295  GaussPenta *gauss1 = NULL;
296  GaussPenta *gauss2 = NULL;
297  gauss1=new GaussPenta(0,1,2,order);
298  gauss2=new GaussPenta(0,1,2,order);
299 
300  xy_list[0][0]=r1; xy_list[0][1]=0;
301  xy_list[1][0]=0; xy_list[1][1]=1.;
302  xy_list[2][0]=0; xy_list[2][1]=r2;
303 
304  //gauss1->Echo();
305  for(ig=0;ig<gauss1->numgauss;ig++){
306  x = gauss1->coords1[ig]*xy_list[0][0] + gauss1->coords2[ig]*xy_list[1][0] + gauss1->coords3[ig]*xy_list[2][0];
307  y = gauss1->coords1[ig]*xy_list[0][1] + gauss1->coords2[ig]*xy_list[1][1] + gauss1->coords3[ig]*xy_list[2][1];
308 
309  switch(index){
310  case 0:
311  gauss1->coords1[ig] = 1.-x-y;
312  gauss1->coords2[ig] = x;
313  gauss1->coords3[ig] = y;
314  break;
315  case 1:
316  gauss1->coords1[ig] = y;
317  gauss1->coords2[ig] = 1.-x-y;
318  gauss1->coords3[ig] = x;
319  break;
320  case 2:
321  gauss1->coords1[ig] = x;
322  gauss1->coords2[ig] = y;
323  gauss1->coords3[ig] = 1.-x-y;
324  break;
325  default:
326  _error_("index "<<index<<" not supported yet");
327  }
328  gauss1->weights[ig] = gauss1->weights[ig]*r1*(1-r2);
329  }
330  //gauss1->Echo();
331  xy_list[0][0]=r1; xy_list[0][1]=0;
332  xy_list[1][0]=1.; xy_list[1][1]=0;
333  xy_list[2][0]=0; xy_list[2][1]=1.;
334 
335  //gauss2->Echo();
336  for(ig=0;ig<gauss2->numgauss;ig++){
337  x = gauss2->coords1[ig]*xy_list[0][0] + gauss2->coords2[ig]*xy_list[1][0] + gauss2->coords3[ig]*xy_list[2][0];
338  y = gauss2->coords1[ig]*xy_list[0][1] + gauss2->coords2[ig]*xy_list[1][1] + gauss2->coords3[ig]*xy_list[2][1];
339 
340  switch(index){
341  case 0:
342  gauss2->coords1[ig] = 1.-x-y;
343  gauss2->coords2[ig] = x;
344  gauss2->coords3[ig] = y;
345  break;
346  case 1:
347  gauss2->coords1[ig] = y;
348  gauss2->coords2[ig] = 1.-x-y;
349  gauss2->coords3[ig] = x;
350  break;
351  case 2:
352  gauss2->coords1[ig] = x;
353  gauss2->coords2[ig] = y;
354  gauss2->coords3[ig] = 1.-x-y;
355  break;
356  default:
357  _error_("index "<<index<<" not supported yet");
358  }
359  gauss2->weights[ig] = gauss2->weights[ig]*(1-r1);
360  }
361 
362  this->numgauss = gauss1->numgauss + gauss2->numgauss;
363  this->coords1=xNew<IssmDouble>(this->numgauss);
364  this->coords2=xNew<IssmDouble>(this->numgauss);
365  this->coords3=xNew<IssmDouble>(this->numgauss);
366  this->coords4=xNew<IssmDouble>(this->numgauss);
367  this->weights=xNew<IssmDouble>(this->numgauss);
368 
369  for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points
370  this->coords1[ig]=gauss1->coords1[ig];
371  this->coords2[ig]=gauss1->coords2[ig];
372  this->coords3[ig]=gauss1->coords3[ig];
373  this->coords4[ig]=gauss1->coords4[ig];
374  this->weights[ig]=gauss1->weights[ig];
375  }
376  for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points
377  this->coords1[gauss1->numgauss+ig]=gauss2->coords1[ig];
378  this->coords2[gauss1->numgauss+ig]=gauss2->coords2[ig];
379  this->coords3[gauss1->numgauss+ig]=gauss2->coords3[ig];
380  this->coords4[gauss1->numgauss+ig]=gauss2->coords4[ig];
381  this->weights[gauss1->numgauss+ig]=gauss2->weights[ig];
382  }
383 
384  /*Delete gauss points*/
385  delete gauss1;
386  delete gauss2;
387  }
388 
389  /*Initialize static fields as undefined*/
390  weight=UNDEF;
391  coord1=UNDEF;
392  coord2=UNDEF;
393  coord3=UNDEF;
394  coord4=UNDEF;
395 }

◆ GaussPenta() [7/8]

GaussPenta::GaussPenta ( IssmDouble  area_coordinates[4][3],
int  order_horiz,
int  order_vert 
)

Definition at line 397 of file GaussPenta.cpp.

397  {/*{{{*/
398 
399  /*Intermediaties*/
400  IssmPDouble *seg_horiz_coords = NULL;
401  IssmPDouble *seg_horiz_weights = NULL;
402  IssmPDouble *seg_vert_coords = NULL;
403  IssmPDouble *seg_vert_weights = NULL;
404 
405  /*get the gauss points using the product of two line rules*/
406  GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
407  GaussLegendreLinear(&seg_vert_coords, &seg_vert_weights, order_vert);
408 
409  /*Allocate GaussPenta fields*/
410  numgauss=order_horiz*order_vert;
411  coords1=xNew<IssmDouble>(numgauss);
412  coords2=xNew<IssmDouble>(numgauss);
413  coords3=xNew<IssmDouble>(numgauss);
414  coords4=xNew<IssmDouble>(numgauss);
415  weights=xNew<IssmDouble>(numgauss);
416 
417  /*Quads: get the gauss points using the product of two line rules */
418  for(int i=0;i<order_horiz;i++){
419  for(int j=0;j<order_vert;j++){
420  coords1[i*order_vert+j]=0.5*(area_coordinates[0][0]+area_coordinates[1][0]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][0]-area_coordinates[0][0]);
421  coords2[i*order_vert+j]=0.5*(area_coordinates[0][1]+area_coordinates[1][1]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][1]-area_coordinates[0][1]);
422  coords3[i*order_vert+j]=0.5*(area_coordinates[0][2]+area_coordinates[1][2]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][2]-area_coordinates[0][2]);
423  coords4[i*order_vert+j]=seg_vert_coords[j];
424  weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
425  }
426  }
427 
428  /*clean-up*/
429  xDelete<IssmPDouble>(seg_horiz_coords);
430  xDelete<IssmPDouble>(seg_horiz_weights);
431  xDelete<IssmPDouble>(seg_vert_coords);
432  xDelete<IssmPDouble>(seg_vert_weights);
433 }

◆ GaussPenta() [8/8]

GaussPenta::GaussPenta ( IssmDouble  area_coordinates[2][3],
int  order_horiz 
)

Definition at line 435 of file GaussPenta.cpp.

435  {/*{{{*/
436 
437  /*Intermediaties*/
438  IssmPDouble *seg_horiz_coords = NULL;
439  IssmPDouble *seg_horiz_weights = NULL;
440 
441  /*get the gauss points using the product of two line rules*/
442  GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
443 
444  /*Allocate GaussPenta fields*/
445  numgauss=order_horiz;
446  coords1=xNew<IssmDouble>(numgauss);
447  coords2=xNew<IssmDouble>(numgauss);
448  coords3=xNew<IssmDouble>(numgauss);
449  coords4=xNew<IssmDouble>(numgauss);
450  weights=xNew<IssmDouble>(numgauss);
451 
452  /*Quads: get the gauss points using the product of two line rules */
453  for(int i=0;i<order_horiz;i++){
454  coords1[i]=0.5*(area_coordinates[0][0]+area_coordinates[1][0]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][0]-area_coordinates[0][0]);
455  coords2[i]=0.5*(area_coordinates[0][1]+area_coordinates[1][1]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][1]-area_coordinates[0][1]);
456  coords3[i]=0.5*(area_coordinates[0][2]+area_coordinates[1][2]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][2]-area_coordinates[0][2]);
457  coords4[i]=0.;
458  weights[i]=seg_horiz_weights[i];
459  }
460 
461  /*clean-up*/
462  xDelete<IssmPDouble>(seg_horiz_coords);
463  xDelete<IssmPDouble>(seg_horiz_weights);
464 }

◆ ~GaussPenta()

GaussPenta::~GaussPenta ( )

Definition at line 466 of file GaussPenta.cpp.

466  {/*{{{*/
467  xDelete<IssmDouble>(weights);
468  xDelete<IssmDouble>(coords1);
469  xDelete<IssmDouble>(coords2);
470  xDelete<IssmDouble>(coords3);
471  xDelete<IssmDouble>(coords4);
472 }

Member Function Documentation

◆ begin()

int GaussPenta::begin ( void  )
virtual

Implements Gauss.

Definition at line 476 of file GaussPenta.cpp.

476  {/*{{{*/
477 
478  /*Check that this has been initialized*/
479  _assert_(numgauss>0);
480  _assert_(weights);
481  _assert_(coords1);
482  _assert_(coords2);
483  _assert_(coords3);
484  _assert_(coords4);
485 
486  /*return first gauss index*/
487  return 0;
488 }

◆ Echo()

void GaussPenta::Echo ( void  )
virtual

Implements Gauss.

Definition at line 490 of file GaussPenta.cpp.

490  {/*{{{*/
491 
492  _printf_("GaussPenta:\n");
493  _printf_(" numgauss: " << numgauss << "\n");
494 
495  if (weights){
496  _printf_(" weights = [");
497  for(int i=0;i<numgauss;i++) _printf_(" " << weights[i] << "\n");
498  _printf_("]\n");
499  }
500  else _printf_("weights = NULL\n");
501  if (coords1){
502  _printf_(" coords1 = [");
503  for(int i=0;i<numgauss;i++) _printf_(" " << coords1[i] << "\n");
504  _printf_("]\n");
505  }
506  else _printf_("coords1 = NULL\n");
507  if (coords2){
508  _printf_(" coords2 = [");
509  for(int i=0;i<numgauss;i++) _printf_(" " << coords2[i] << "\n");
510  _printf_("]\n");
511  }
512  else _printf_("coords2 = NULL\n");
513  if (coords3){
514  _printf_(" coords3 = [");
515  for(int i=0;i<numgauss;i++) _printf_(" " << coords3[i] << "\n");
516  _printf_("]\n");
517  }
518  else _printf_("coords3 = NULL\n");
519  if (coords4){
520  _printf_(" coords4 = [");
521  for(int i=0;i<numgauss;i++) _printf_(" " << coords4[i] << "\n");
522  _printf_("]\n");
523  }
524  else _printf_("coords4 = NULL\n");
525 
526  _printf_(" weight = " << weight << "\n");
527  _printf_(" coord1 = " << coord1 << "\n");
528  _printf_(" coord2 = " << coord2 << "\n");
529  _printf_(" coord3 = " << coord3 << "\n");
530  _printf_(" coord4 = " << coord4 << "\n");
531 
532 }

◆ end()

int GaussPenta::end ( void  )
virtual

Implements Gauss.

Definition at line 534 of file GaussPenta.cpp.

534  {/*{{{*/
535 
536  /*Check that this has been initialized*/
537  _assert_(numgauss>0);
538  _assert_(weights);
539  _assert_(coords1);
540  _assert_(coords2);
541  _assert_(coords3);
542  _assert_(coords4);
543 
544  /*return last gauss index +1*/
545  return numgauss;
546 }

◆ Enum()

int GaussPenta::Enum ( void  )
virtual

Implements Gauss.

Definition at line 548 of file GaussPenta.cpp.

548  {/*{{{*/
549  return GaussPentaEnum;
550 }

◆ GaussFaceTria()

void GaussPenta::GaussFaceTria ( int  index1,
int  index2,
int  index3,
int  order 
)

Definition at line 552 of file GaussPenta.cpp.

552  {/*{{{*/
553 
554  /*in debugging mode: check that the default constructor has been called*/
555  _assert_(numgauss==-1);
556 
557  /*Basal Tria*/
558  if(index1==0 && index2==1 && index3==2){
559  GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
560  coords4=xNew<IssmDouble>(numgauss);
561  for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
562  }
563  else{
564  _error_("Tria not supported yet");
565  }
566 
567 }

◆ GaussNode()

void GaussPenta::GaussNode ( int  finitelement,
int  iv 
)
virtual

Implements Gauss.

Definition at line 583 of file GaussPenta.cpp.

583  {/*{{{*/
584 
585  /*in debugging mode: check that the default constructor has been called*/
586  _assert_(numgauss==-1);
587 
588  /*update static arrays*/
589  switch(finiteelement){
590  case P1Enum: case P1DGEnum:
591  switch(iv){
592  case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
593  case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
594  case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
595  case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
596  case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
597  case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
598  default: _error_("node index should be in [0 5]");
599  }
600  break;
601  case P1xP2Enum:
602  switch(iv){
603  case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
604  case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
605  case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
606  case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
607  case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
608  case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
609 
610  case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
611  case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
612  case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
613  default: _error_("node index should be in [0 8]");
614  }
615  break;
616  case P1xP3Enum:
617  switch(iv){
618  case 0 : coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
619  case 1 : coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
620  case 2 : coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
621  case 3 : coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
622  case 4 : coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
623  case 5 : coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
624 
625  case 6 : coord1=1.; coord2=0.; coord3=0.; coord4=-1./3.; break;
626  case 7 : coord1=0.; coord2=1.; coord3=0.; coord4=-1./3.; break;
627  case 8 : coord1=0.; coord2=0.; coord3=1.; coord4=-1./3.; break;
628  case 9 : coord1=1.; coord2=0.; coord3=0.; coord4=+1./3.; break;
629  case 10: coord1=0.; coord2=1.; coord3=0.; coord4=+1./3.; break;
630  case 11: coord1=0.; coord2=0.; coord3=1.; coord4=+1./3.; break;
631  default: _error_("node index should be in [0 11]");
632  }
633  break;
634  case P1xP4Enum:
635  switch(iv){
636  case 0 : coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
637  case 1 : coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
638  case 2 : coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
639  case 3 : coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
640  case 4 : coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
641  case 5 : coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
642 
643  case 6 : coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
644  case 7 : coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
645  case 8 : coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
646 
647  case 9 : coord1=1.; coord2=0.; coord3=0.; coord4=-0.5; break;
648  case 10: coord1=0.; coord2=1.; coord3=0.; coord4=-0.5; break;
649  case 11: coord1=0.; coord2=0.; coord3=1.; coord4=-0.5; break;
650 
651  case 12: coord1=1.; coord2=0.; coord3=0.; coord4=+0.5; break;
652  case 13: coord1=0.; coord2=1.; coord3=0.; coord4=+0.5; break;
653  case 14: coord1=0.; coord2=0.; coord3=1.; coord4=+0.5; break;
654  default: _error_("node index should be in [0 14]");
655  }
656  break;
657  case P2xP1Enum:
658  switch(iv){
659  case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
660  case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
661  case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
662  case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
663  case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
664  case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
665 
666  case 6: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
667  case 7: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
668  case 8: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
669  case 9: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
670  case 10: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
671  case 11: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
672  default: _error_("node index should be in [0 11]");
673  }
674  break;
675  case P1bubbleEnum: case P1bubblecondensedEnum:
676  switch(iv){
677  case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
678  case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
679  case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
680  case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
681  case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
682  case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
683  case 6: coord1=1./3.; coord2=1./3.; coord3=1./3.; coord4=0.; break;
684  default: _error_("node index should be in [0 6]");
685  }
686  break;
687  case P2Enum:
688  switch(iv){
689  case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
690  case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
691  case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
692  case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
693  case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
694  case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
695 
696  case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
697  case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
698  case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
699 
700  case 9: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
701  case 10: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
702  case 11: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
703  case 12: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
704  case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
705  case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
706 
707  case 15: coord1=0.; coord2=.5; coord3=.5; coord4=0.;break;
708  case 16: coord1=.5; coord2=0.; coord3=.5; coord4=0.;break;
709  case 17: coord1=.5; coord2=.5; coord3=0.; coord4=0.;break;
710  default: _error_("node index should be in [0 17]");
711  }
712  break;
713  case P2bubbleEnum: case P2bubblecondensedEnum:
714  switch(iv){
715  case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
716  case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
717  case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
718  case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
719  case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
720  case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
721 
722  case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
723  case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
724  case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
725 
726  case 9: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
727  case 10: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
728  case 11: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
729  case 12: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
730  case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
731  case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
732 
733  case 15: coord1=0.; coord2=.5; coord3=.5; coord4=0.;break;
734  case 16: coord1=.5; coord2=0.; coord3=.5; coord4=0.;break;
735  case 17: coord1=.5; coord2=.5; coord3=0.; coord4=0.;break;
736 
737  case 18: coord1=1./3.; coord2=1./3.; coord3=1./3.; coord4=0.; break;
738  default: _error_("node index should be in [0 18]");
739  }
740  break;
741  case P2xP4Enum:
742  switch(iv){
743  case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
744  case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
745  case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
746  case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
747  case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
748  case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
749 
750  case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
751  case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
752  case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
753 
754  case 9: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
755  case 10: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
756  case 11: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
757  case 12: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
758  case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
759  case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
760 
761  case 15: coord1=1.; coord2=0.; coord3=0.; coord4=-.5; break;
762  case 16: coord1=0.; coord2=1.; coord3=0.; coord4=-.5; break;
763  case 17: coord1=0.; coord2=0.; coord3=1.; coord4=-.5; break;
764  case 18: coord1=1.; coord2=0.; coord3=0.; coord4=+.5; break;
765  case 19: coord1=0.; coord2=1.; coord3=0.; coord4=+.5; break;
766  case 20: coord1=0.; coord2=0.; coord3=1.; coord4=+.5; break;
767 
768  case 21: coord1=0.; coord2=.5; coord3=.5; coord4=-.5;break;
769  case 22: coord1=.5; coord2=0.; coord3=.5; coord4=-.5;break;
770  case 23: coord1=.5; coord2=.5; coord3=0.; coord4=-.5;break;
771  case 24: coord1=0.; coord2=.5; coord3=.5; coord4=0.;break;
772  case 25: coord1=.5; coord2=0.; coord3=.5; coord4=0.;break;
773  case 26: coord1=.5; coord2=.5; coord3=0.; coord4=0.;break;
774  case 27: coord1=0.; coord2=.5; coord3=.5; coord4=+.5;break;
775  case 28: coord1=.5; coord2=0.; coord3=.5; coord4=+.5;break;
776  case 29: coord1=.5; coord2=.5; coord3=0.; coord4=+.5;break;
777  default: _error_("node index should be in [0 29]");
778  }
779  break;
780  default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
781  }
782 
783 }

◆ GaussPoint()

void GaussPenta::GaussPoint ( int  ig)
virtual

Implements Gauss.

Definition at line 569 of file GaussPenta.cpp.

569  {/*{{{*/
570 
571  /*Check input in debugging mode*/
572  _assert_(ig>=0 && ig< numgauss);
573 
574  /*update static arrays*/
575  weight=weights[ig];
576  coord1=coords1[ig];
577  coord2=coords2[ig];
578  coord3=coords3[ig];
579  coord4=coords4[ig];
580 
581 }

◆ GaussVertex()

void GaussPenta::GaussVertex ( int  iv)
virtual

Implements Gauss.

Definition at line 785 of file GaussPenta.cpp.

785  {/*{{{*/
786 
787  /*in debugging mode: check that the default constructor has been called*/
788  _assert_(numgauss==-1);
789 
790  /*update static arrays*/
791  switch(iv){
792  case 0: coord1=1.; coord2=0.; coord3=0.; coord4= -1.; break;
793  case 1: coord1=0.; coord2=1.; coord3=0.; coord4= -1.; break;
794  case 2: coord1=0.; coord2=0.; coord3=1.; coord4= -1.; break;
795  case 3: coord1=1.; coord2=0.; coord3=0.; coord4= +1.; break;
796  case 4: coord1=0.; coord2=1.; coord3=0.; coord4= +1.; break;
797  case 5: coord1=0.; coord2=0.; coord3=1.; coord4= +1.; break;
798  default: _error_("vertex index should be in [0 5]");
799 
800  }
801 
802 }

◆ SynchronizeGaussBase()

void GaussPenta::SynchronizeGaussBase ( Gauss *  gauss)
virtual

Implements Gauss.

Definition at line 804 of file GaussPenta.cpp.

804  {/*{{{*/
805 
806  _assert_(gauss->Enum()==GaussTriaEnum);
807  GaussTria* gauss_tria = xDynamicCast<GaussTria*>(gauss);
808 
809  gauss_tria->coord1=this->coord1;
810  gauss_tria->coord2=this->coord2;
811  gauss_tria->coord3=this->coord3;
812  gauss_tria->weight=UNDEF;
813 }

Field Documentation

◆ numgauss

int GaussPenta::numgauss
private

Definition at line 16 of file GaussPenta.h.

◆ weights

IssmDouble* GaussPenta::weights
private

Definition at line 17 of file GaussPenta.h.

◆ coords1

IssmDouble* GaussPenta::coords1
private

Definition at line 18 of file GaussPenta.h.

◆ coords2

IssmDouble* GaussPenta::coords2
private

Definition at line 19 of file GaussPenta.h.

◆ coords3

IssmDouble* GaussPenta::coords3
private

Definition at line 20 of file GaussPenta.h.

◆ coords4

IssmDouble* GaussPenta::coords4
private

Definition at line 21 of file GaussPenta.h.

◆ coord1

IssmDouble GaussPenta::coord1

Definition at line 24 of file GaussPenta.h.

◆ coord2

IssmDouble GaussPenta::coord2

Definition at line 25 of file GaussPenta.h.

◆ coord3

IssmDouble GaussPenta::coord3

Definition at line 26 of file GaussPenta.h.

◆ coord4

IssmDouble GaussPenta::coord4

Definition at line 27 of file GaussPenta.h.

The documentation for this class was generated from the following files:
GaussPenta::GaussPenta
GaussPenta()
Definition: GaussPenta.cpp:16
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
GaussPenta::coord1
IssmDouble coord1
Definition: GaussPenta.h:24
GaussPenta::coords2
IssmDouble * coords2
Definition: GaussPenta.h:19
GaussTria::coord1
IssmDouble coord1
Definition: GaussTria.h:22
GaussPenta::coord3
IssmDouble coord3
Definition: GaussPenta.h:26
P1DGEnum
@ P1DGEnum
Definition: EnumDefinitions.h:1215
P2xP1Enum
@ P2xP1Enum
Definition: EnumDefinitions.h:1226
P1bubblecondensedEnum
@ P1bubblecondensedEnum
Definition: EnumDefinitions.h:1219
GaussTria::coord3
IssmDouble coord3
Definition: GaussTria.h:24
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
GaussTria
Definition: GaussTria.h:12
GaussTria::coord2
IssmDouble coord2
Definition: GaussTria.h:23
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
GaussPenta
Definition: GaussPenta.h:13
GaussPenta::coord2
IssmDouble coord2
Definition: GaussPenta.h:25
GaussPenta::weights
IssmDouble * weights
Definition: GaussPenta.h:17
UNDEF
#define UNDEF
Definition: constants.h:8
P2bubbleEnum
@ P2bubbleEnum
Definition: EnumDefinitions.h:1224
P2xP4Enum
@ P2xP4Enum
Definition: EnumDefinitions.h:1227
GaussLegendreLinear
void GaussLegendreLinear(IssmPDouble **pxgaus, IssmPDouble **pxwgt, int ngaus)
Definition: GaussPoints.cpp:11
Gauss::Enum
virtual int Enum(void)=0
P1bubbleEnum
@ P1bubbleEnum
Definition: EnumDefinitions.h:1218
GaussPenta::numgauss
int numgauss
Definition: GaussPenta.h:16
GaussPenta::coords4
IssmDouble * coords4
Definition: GaussPenta.h:21
GaussPenta::coord4
IssmDouble coord4
Definition: GaussPenta.h:27
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
P2bubblecondensedEnum
@ P2bubblecondensedEnum
Definition: EnumDefinitions.h:1225
GaussPenta::coords1
IssmDouble * coords1
Definition: GaussPenta.h:18
GaussTriaEnum
@ GaussTriaEnum
Definition: EnumDefinitions.h:1081
GaussLegendreTria
void GaussLegendreTria(int *pngaus, IssmDouble **pl1, IssmDouble **pl2, IssmDouble **pl3, IssmDouble **pwgt, int iord)
Definition: GaussPoints.cpp:95
P2Enum
@ P2Enum
Definition: EnumDefinitions.h:1223
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
P1xP4Enum
@ P1xP4Enum
Definition: EnumDefinitions.h:1222
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
P1xP3Enum
@ P1xP3Enum
Definition: EnumDefinitions.h:1221
P1xP2Enum
@ P1xP2Enum
Definition: EnumDefinitions.h:1220
GaussPenta::coords3
IssmDouble * coords3
Definition: GaussPenta.h:20
GaussPentaEnum
@ GaussPentaEnum
Definition: EnumDefinitions.h:1078
Generated on Thu Jul 2 2020 08:09:27 for Ice Sheet System Model by   doxygen 1.8.19