Ice Sheet System Model  4.18
Code documentation
GaussPenta.cpp
Go to the documentation of this file.
1 
5 #include "./GaussPenta.h"
6 #include "./GaussTria.h"
7 #include "../../shared/io/Print/Print.h"
8 #include "../../shared/Exceptions/exceptions.h"
9 #include "../../shared/MemOps/MemOps.h"
10 #include "../../shared/Numerics/recast.h"
11 #include "../../shared/Enum/Enum.h"
12 #include "../../shared/Numerics/GaussPoints.h"
13 #include "../../shared/Numerics/constants.h"
14 
15 /*GaussPenta constructors and destructors:*/
16 GaussPenta::GaussPenta(){/*{{{*/
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 }
32 /*}}}*/
33 GaussPenta::GaussPenta(int order_horiz,int order_vert){/*{{{*/
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 }
85 /*}}}*/
86 GaussPenta::GaussPenta(int index1, int index2,int order){/*{{{*/
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 }
141 /*}}}*/
142 GaussPenta::GaussPenta(int index1, int index2, int index3, int order){/*{{{*/
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 }
169 /*}}}*/
170 GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){/*{{{*/
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 }
235 /*}}}*/
236 GaussPenta::GaussPenta(int index,IssmDouble r1,IssmDouble r2,bool mainlyfloating,int order){/*{{{*/
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 }
396 /*}}}*/
397 GaussPenta::GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert){/*{{{*/
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 }
434 /*}}}*/
435 GaussPenta::GaussPenta(IssmDouble area_coordinates[2][3],int order_horiz){/*{{{*/
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 }
465 /*}}}*/
466 GaussPenta::~GaussPenta(){/*{{{*/
467  xDelete<IssmDouble>(weights);
468  xDelete<IssmDouble>(coords1);
469  xDelete<IssmDouble>(coords2);
470  xDelete<IssmDouble>(coords3);
471  xDelete<IssmDouble>(coords4);
472 }
473 /*}}}*/
474 
475 /*Methods*/
476 int GaussPenta::begin(void){/*{{{*/
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 }
489 /*}}}*/
490 void GaussPenta::Echo(void){/*{{{*/
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 }
533 /*}}}*/
534 int GaussPenta::end(void){/*{{{*/
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 }
547 /*}}}*/
548 int GaussPenta::Enum(void){/*{{{*/
549  return GaussPentaEnum;
550 }
551 /*}}}*/
552 void GaussPenta::GaussFaceTria(int index1, int index2, int index3, int order){/*{{{*/
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 }
568 /*}}}*/
569 void GaussPenta::GaussPoint(int ig){/*{{{*/
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 }
582 /*}}}*/
583 void GaussPenta::GaussNode(int finiteelement,int iv){/*{{{*/
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 }
784 /*}}}*/
785 void GaussPenta::GaussVertex(int iv){/*{{{*/
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 }
803 /*}}}*/
804 void GaussPenta::SynchronizeGaussBase(Gauss* gauss){/*{{{*/
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 }
814 /*}}}*/
GaussPenta::Enum
int Enum(void)
Definition: GaussPenta.cpp:548
GaussPenta::GaussPenta
GaussPenta()
Definition: GaussPenta.cpp:16
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
GaussPenta::GaussPoint
void GaussPoint(int ig)
Definition: GaussPenta.cpp:569
GaussPenta::GaussVertex
void GaussVertex(int iv)
Definition: GaussPenta.cpp:785
GaussPenta.h
: header file for node object
IssmDouble
double IssmDouble
Definition: types.h:37
GaussPenta::SynchronizeGaussBase
void SynchronizeGaussBase(Gauss *gauss)
Definition: GaussPenta.cpp:804
_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
GaussPenta::Echo
void Echo(void)
Definition: GaussPenta.cpp:490
P2xP1Enum
@ P2xP1Enum
Definition: EnumDefinitions.h:1226
GaussPenta::~GaussPenta
~GaussPenta()
Definition: GaussPenta.cpp:466
P1bubblecondensedEnum
@ P1bubblecondensedEnum
Definition: EnumDefinitions.h:1219
GaussTria::coord3
IssmDouble coord3
Definition: GaussTria.h:24
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
GaussPenta::begin
int begin(void)
Definition: GaussPenta.cpp:476
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
GaussTria.h
: header file for node object
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
GaussPenta::end
int end(void)
Definition: GaussPenta.cpp:534
GaussPenta::GaussFaceTria
void GaussFaceTria(int index1, int index2, int index3, int order)
Definition: GaussPenta.cpp:552
P2bubblecondensedEnum
@ P2bubblecondensedEnum
Definition: EnumDefinitions.h:1225
GaussPenta::GaussNode
void GaussNode(int finitelement, int iv)
Definition: GaussPenta.cpp:583
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
Gauss
Definition: Gauss.h:8
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:17 for Ice Sheet System Model by   doxygen 1.8.19