source: issm/trunk/src/c/objects/Gauss/GaussPenta.cpp@ 5808

Last change on this file since 5808 was 5808, checked in by Mathieu Morlighem, 15 years ago

New Icefront MUCH cleaner. Enjoy
File size: 10.5 KB
Line 
1/*!\file GaussPenta.c
2 * \brief: implementation of the GaussPenta object
3 */
4
5/*Include files: {{{1*/
6#include "./../objects.h"
7/*}}}*/
8
9/*GaussPenta constructors and destructors:*/
10/*FUNCTION GaussPenta::GaussPenta() {{{1*/
11GaussPenta::GaussPenta(){
12
13 numgauss=-1;
14
15 weights=NULL;
16 coords1=NULL;
17 coords2=NULL;
18 coords3=NULL;
19 coords4=NULL;
20
21 weight=UNDEF;
22 coord1=UNDEF;
23 coord2=UNDEF;
24 coord3=UNDEF;
25 coord4=UNDEF;
26}
27/*}}}*/
28/*FUNCTION GaussPenta::GaussPenta(int order_horiz,int order_vert) {{{1*/
29GaussPenta::GaussPenta(int order_horiz,int order_vert){
30
31 /*Intermediaries*/
32 int ighoriz,igvert;
33 int numgauss_horiz;
34 int numgauss_vert;
35 double *coords1_horiz = NULL;
36 double *coords2_horiz = NULL;
37 double *coords3_horiz = NULL;
38 double *weights_horiz = NULL;
39 double *coords_vert = NULL;
40 double *weights_vert = NULL;
41
42 /*Get gauss points*/
43 GaussLegendreTria(&numgauss_horiz,&coords1_horiz,&coords2_horiz,&coords3_horiz,&weights_horiz,order_horiz);
44 GaussLegendreLinear(&coords_vert,&weights_vert,order_vert);
45 numgauss_vert=order_vert;
46
47 /*Allocate GaussPenta fields*/
48 numgauss=numgauss_horiz*numgauss_vert;
49 coords1=(double*)xmalloc(numgauss*sizeof(double));
50 coords2=(double*)xmalloc(numgauss*sizeof(double));
51 coords3=(double*)xmalloc(numgauss*sizeof(double));
52 coords4=(double*)xmalloc(numgauss*sizeof(double));
53 weights=(double*)xmalloc(numgauss*sizeof(double));
54
55 /*Combine Horizontal and vertical points*/
56 for (ighoriz=0; ighoriz<numgauss_horiz; ighoriz++){
57 for (igvert=0; igvert<numgauss_vert; igvert++){
58 coords1[numgauss_vert*ighoriz+igvert]=coords1_horiz[ighoriz];
59 coords2[numgauss_vert*ighoriz+igvert]=coords2_horiz[ighoriz];
60 coords3[numgauss_vert*ighoriz+igvert]=coords3_horiz[ighoriz];
61 coords4[numgauss_vert*ighoriz+igvert]=coords_vert[igvert];
62 weights[numgauss_vert*ighoriz+igvert]=weights_horiz[ighoriz]*weights_vert[igvert];
63 }
64 }
65
66 /*Initialize static fields as undefinite*/
67 weight=UNDEF;
68 coord1=UNDEF;
69 coord2=UNDEF;
70 coord3=UNDEF;
71 coord4=UNDEF;
72
73 /*Clean up*/
74 xfree((void**)&coords1_horiz);
75 xfree((void**)&coords2_horiz);
76 xfree((void**)&coords3_horiz);
77 xfree((void**)&coords_vert);
78 xfree((void**)&weights_horiz);
79 xfree((void**)&weights_vert);
80}
81/*}}}*/
82/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int order){{{1*/
83GaussPenta::GaussPenta(int index1, int index2,int order){
84
85 /*Intermediaties*/
86 double *seg_coords = NULL;
87 double *seg_weights = NULL;
88 int i;
89
90 /*Get Segment gauss points*/
91 numgauss=order;
92 GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);
93
94 /*Allocate GaussPenta fields*/
95 coords1=(double*)xmalloc(numgauss*sizeof(double));
96 coords2=(double*)xmalloc(numgauss*sizeof(double));
97 coords3=(double*)xmalloc(numgauss*sizeof(double));
98 coords4=(double*)xmalloc(numgauss*sizeof(double));
99 weights=(double*)xmalloc(numgauss*sizeof(double));
100
101 if(index1==0 && index2==3){
102 for(i=0;i<numgauss;i++) coords1[i]=1.0;
103 for(i=0;i<numgauss;i++) coords2[i]=0.0;
104 for(i=0;i<numgauss;i++) coords3[i]=0.0;
105 for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
106 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
107 }
108 else if (index1==1 && index2==4){
109 for(i=0;i<numgauss;i++) coords1[i]=0.0;
110 for(i=0;i<numgauss;i++) coords2[i]=1.0;
111 for(i=0;i<numgauss;i++) coords3[i]=0.0;
112 for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
113 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
114 }
115 else if (index1==2 && index2==5){
116 for(i=0;i<numgauss;i++) coords1[i]=0.0;
117 for(i=0;i<numgauss;i++) coords2[i]=0.0;
118 for(i=0;i<numgauss;i++) coords3[i]=1.0;
119 for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
120 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
121 }
122 else{
123 ISSMERROR("Penta not supported yet");
124 }
125
126 /*Initialize static fields as undefined*/
127 weight=UNDEF;
128 coord1=UNDEF;
129 coord2=UNDEF;
130 coord3=UNDEF;
131 coord4=UNDEF;
132
133 /*clean up*/
134 xfree((void**)&seg_coords);
135 xfree((void**)&seg_weights);
136
137}
138/*}}}*/
139/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int index3, int order){{{1*/
140GaussPenta::GaussPenta(int index1, int index2, int index3, int order){
141
142 /*Basal Tria*/
143 if(index1==0 && index2==1 && index3==2){
144
145 /*Get GaussTria*/
146 GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
147
148 /*compute z coordinate*/
149 coords4=(double*)xmalloc(numgauss*sizeof(double));
150 for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
151 }
152 else{
153 ISSMERROR("Tria not supported yet");
154 }
155
156}
157/*}}}*/
158/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){{{1*/
159GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){
160
161 /*Intermediaties*/
162 double *seg_horiz_coords = NULL;
163 double *seg_horiz_weights = NULL;
164 double *seg_vert_coords = NULL;
165 double *seg_vert_weights = NULL;
166 int i,j;
167
168 /*get the gauss points using the product of two line rules*/
169 GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
170 GaussLegendreLinear(&seg_vert_coords, &seg_vert_weights, order_vert);
171
172 /*Allocate GaussPenta fields*/
173 numgauss=order_horiz*order_vert;
174 coords1=(double*)xmalloc(numgauss*sizeof(double));
175 coords2=(double*)xmalloc(numgauss*sizeof(double));
176 coords3=(double*)xmalloc(numgauss*sizeof(double));
177 coords4=(double*)xmalloc(numgauss*sizeof(double));
178 weights=(double*)xmalloc(numgauss*sizeof(double));
179
180 /*Quads: get the gauss points using the product of two line rules */
181 if(index1==0 && index2==1 && index3==4 && index4==3){
182 for(i=0;i<order_horiz;i++){
183 for(j=0;j<order_vert;j++){
184 coords1[i*order_vert+j]= 0.5*(1-seg_horiz_coords[i]);
185 coords2[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
186 coords3[i*order_vert+j]=0.0;
187 coords4[i*order_vert+j]=seg_vert_coords[j];
188 weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
189 }
190 }
191 }
192 else if(index1==1 && index2==2 && index3==5 && index4==4){
193 for(i=0;i<order_horiz;i++){
194 for(j=0;j<order_vert;j++){
195 coords1[i*order_vert+j]=0.0;
196 coords2[i*order_vert+j]= 0.5*(1-seg_horiz_coords[i]);
197 coords3[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
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==2 && index2==0 && index3==3 && index4==5){
204 for(i=0;i<order_horiz;i++){
205 for(j=0;j<order_vert;j++){
206 coords1[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
207 coords2[i*order_vert+j]=0.0;
208 coords3[i*order_vert+j]= 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{
215 ISSMERROR("Tria not supported yet (user provided indices %i %i %i %i)",index1,index2,index3,index4);
216 }
217
218 /*clean-up*/
219 xfree((void**)&seg_horiz_coords);
220 xfree((void**)&seg_horiz_weights);
221 xfree((void**)&seg_vert_coords);
222 xfree((void**)&seg_vert_weights);
223}
224/*}}}*/
225/*FUNCTION GaussPenta::~GaussPenta(){{{1*/
226GaussPenta::~GaussPenta(){
227 xfree((void**)&weights);
228 xfree((void**)&coords1);
229 xfree((void**)&coords2);
230 xfree((void**)&coords3);
231 xfree((void**)&coords4);
232}
233/*}}}*/
234
235/*Methods*/
236/*FUNCTION GaussPenta::Echo{{{1*/
237void GaussPenta::Echo(void){
238
239 printf("GaussPenta:\n");
240 printf(" numgauss: %i\n",numgauss);
241
242 if (weights){
243 printf(" weights = [");
244 for(int i=0;i<numgauss;i++) printf(" %g\n",weights[i]);
245 printf("]\n");
246 }
247 else printf("weights = NULL\n");
248 if (coords1){
249 printf(" coords1 = [");
250 for(int i=0;i<numgauss;i++) printf(" %g\n",coords1[i]);
251 printf("]\n");
252 }
253 else printf("coords1 = NULL\n");
254 if (coords2){
255 printf(" coords2 = [");
256 for(int i=0;i<numgauss;i++) printf(" %g\n",coords2[i]);
257 printf("]\n");
258 }
259 else printf("coords2 = NULL\n");
260 if (coords3){
261 printf(" coords3 = [");
262 for(int i=0;i<numgauss;i++) printf(" %g\n",coords3[i]);
263 printf("]\n");
264 }
265 else printf("coords3 = NULL\n");
266 if (coords4){
267 printf(" coords4 = [");
268 for(int i=0;i<numgauss;i++) printf(" %g\n",coords4[i]);
269 printf("]\n");
270 }
271 else printf("coords4 = NULL\n");
272
273 printf(" weight = %g\n",weight);
274 printf(" coord1 = %g\n",coord1);
275 printf(" coord2 = %g\n",coord2);
276 printf(" coord3 = %g\n",coord3);
277 printf(" coord4 = %g\n",coord4);
278
279}
280/*}}}*/
281/*FUNCTION GaussPenta::GaussCenter{{{1*/
282void GaussPenta::GaussCenter(void){
283
284 /*update static arrays*/
285 coord1=ONETHIRD;
286 coord2=ONETHIRD;
287 coord3=ONETHIRD;
288 coord4=0.0;
289
290}
291/*}}}*/
292/*FUNCTION GaussPenta::GaussPoint{{{1*/
293void GaussPenta::GaussPoint(int ig){
294
295 /*Check input in debugging mode*/
296 ISSMASSERT(ig>=0 && ig< numgauss);
297
298 /*update static arrays*/
299 weight=weights[ig];
300 coord1=coords1[ig];
301 coord2=coords2[ig];
302 coord3=coords3[ig];
303 coord4=coords4[ig];
304
305}
306/*}}}*/
307/*FUNCTION GaussPenta::GaussVertex{{{1*/
308void GaussPenta::GaussVertex(int iv){
309
310 /*in debugging mode: check that the default constructor has been called*/
311 ISSMASSERT(numgauss==-1);
312
313 /*update static arrays*/
314 switch(iv){
315 case 0:
316 coord1=1; coord2=0; coord3=0; coord4= -1;
317 break;
318 case 1:
319 coord1=0; coord2=1; coord3=0; coord4= -1;
320 break;
321 case 2:
322 coord1=0; coord2=0; coord3=1; coord4= -1;
323 break;
324 case 3:
325 coord1=1; coord2=0; coord3=0; coord4= +1;
326 break;
327 case 4:
328 coord1=0; coord2=1; coord3=0; coord4= +1;
329 break;
330 case 5:
331 coord1=0; coord2=0; coord3=1; coord4= +1;
332 break;
333 default:
334 ISSMERROR("vertex index should be in [0 5]");
335
336 }
337
338}
339/*}}}*/
340/*FUNCTION GaussPenta::GaussFaceTria{{{1*/
341void GaussPenta::GaussFaceTria(int index1, int index2, int index3, int order){
342
343 /*in debugging mode: check that the default constructor has been called*/
344 ISSMASSERT(numgauss==-1);
345
346 /*Basal Tria*/
347 if(index1==0 && index2==1 && index3==2){
348 GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
349 coords4=(double*)xmalloc(numgauss*sizeof(double));
350 for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
351 }
352 else{
353 ISSMERROR("Tria not supported yet");
354 }
355
356}
357/*}}}*/
358/*FUNCTION GaussPenta::begin{{{1*/
359int GaussPenta::begin(void){
360
361 /*Check that this has been initialized*/
362 ISSMASSERT(numgauss>0);
363 ISSMASSERT(weights);
364 ISSMASSERT(coords1);
365 ISSMASSERT(coords2);
366 ISSMASSERT(coords3);
367 ISSMASSERT(coords4);
368
369 /*return first gauss index*/
370 return 0;
371}
372/*}}}*/
373/*FUNCTION GaussPenta::end{{{1*/
374int GaussPenta::end(void){
375
376 /*Check that this has been initialized*/
377 ISSMASSERT(numgauss>0);
378 ISSMASSERT(weights);
379 ISSMASSERT(coords1);
380 ISSMASSERT(coords2);
381 ISSMASSERT(coords3);
382 ISSMASSERT(coords4);
383
384 /*return last gauss index +1*/
385 return numgauss;
386}
387/*}}}*/
388/*FUNCTION GaussPenta::SynchronizeGaussTria{{{1*/
389void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){
390
391 gauss_tria->coord1=this->coord1;
392 gauss_tria->coord2=this->coord2;
393 gauss_tria->coord3=this->coord3;
394 gauss_tria->weight=UNDEF;
395}
396/*}}}*/
Note: See TracBrowser for help on using the repository browser.