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

Last change on this file since 12706 was 12706, checked in by Mathieu Morlighem, 13 years ago

merged trunk-jpl and trunk for revision 12703
File size: 10.5 KB
RevLine 
[5641]1/*!\file GaussPenta.c
2 * \brief: implementation of the GaussPenta object
3 */
4
[12706]5/*Include files: {{{*/
[5641]6#include "./../objects.h"
7/*}}}*/
8
9/*GaussPenta constructors and destructors:*/
[12706]10/*FUNCTION GaussPenta::GaussPenta() {{{*/
[5641]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/*}}}*/
[12706]28/*FUNCTION GaussPenta::GaussPenta(int order_horiz,int order_vert) {{{*/
[5641]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;
[12706]49 coords1=xNew<double>(numgauss);
50 coords2=xNew<double>(numgauss);
51 coords3=xNew<double>(numgauss);
52 coords4=xNew<double>(numgauss);
53 weights=xNew<double>(numgauss);
[5641]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*/
[12706]74 xDelete<double>(coords1_horiz);
75 xDelete<double>(coords2_horiz);
76 xDelete<double>(coords3_horiz);
77 xDelete<double>(coords_vert);
78 xDelete<double>(weights_horiz);
79 xDelete<double>(weights_vert);
[5641]80}
81/*}}}*/
[12706]82/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int order){{{*/
[5724]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
[5797]94 /*Allocate GaussPenta fields*/
[12706]95 coords1=xNew<double>(numgauss);
96 coords2=xNew<double>(numgauss);
97 coords3=xNew<double>(numgauss);
98 coords4=xNew<double>(numgauss);
99 weights=xNew<double>(numgauss);
[5724]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{
[12706]123 _error2_("Penta not supported yet");
[5724]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*/
[12706]134 xDelete<double>(seg_coords);
135 xDelete<double>(seg_weights);
[5724]136
137}
138/*}}}*/
[12706]139/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int index3, int order){{{*/
[5722]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*/
[12706]149 coords4=xNew<double>(numgauss);
[5722]150 for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
151 }
[8416]152 /*Upper surface Tria*/
153 else if(index1==3 && index2==4 && index3==5){
154
155 /*Get GaussTria*/
156 GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
157
158 /*compute z coordinate*/
[12706]159 coords4=xNew<double>(numgauss);
[8416]160 for(int i=0;i<numgauss;i++) coords4[i]=1.0;
161 }
[5722]162 else{
[12706]163 _error2_("Tria not supported yet");
[5722]164 }
165
166}
167/*}}}*/
[12706]168/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){{{*/
[5797]169GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){
170
171 /*Intermediaties*/
172 double *seg_horiz_coords = NULL;
173 double *seg_horiz_weights = NULL;
174 double *seg_vert_coords = NULL;
175 double *seg_vert_weights = NULL;
176 int i,j;
177
178 /*get the gauss points using the product of two line rules*/
179 GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
180 GaussLegendreLinear(&seg_vert_coords, &seg_vert_weights, order_vert);
181
182 /*Allocate GaussPenta fields*/
183 numgauss=order_horiz*order_vert;
[12706]184 coords1=xNew<double>(numgauss);
185 coords2=xNew<double>(numgauss);
186 coords3=xNew<double>(numgauss);
187 coords4=xNew<double>(numgauss);
188 weights=xNew<double>(numgauss);
[5797]189
190 /*Quads: get the gauss points using the product of two line rules */
[5808]191 if(index1==0 && index2==1 && index3==4 && index4==3){
[5797]192 for(i=0;i<order_horiz;i++){
193 for(j=0;j<order_vert;j++){
[5808]194 coords1[i*order_vert+j]= 0.5*(1-seg_horiz_coords[i]);
195 coords2[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
196 coords3[i*order_vert+j]=0.0;
197 coords4[i*order_vert+j]=seg_vert_coords[j];
198 weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
[5797]199 }
200 }
201 }
[5808]202 else if(index1==1 && index2==2 && index3==5 && index4==4){
[5797]203 for(i=0;i<order_horiz;i++){
204 for(j=0;j<order_vert;j++){
[5808]205 coords1[i*order_vert+j]=0.0;
206 coords2[i*order_vert+j]= 0.5*(1-seg_horiz_coords[i]);
207 coords3[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
208 coords4[i*order_vert+j]=seg_vert_coords[j];
209 weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
[5797]210 }
211 }
212 }
[5808]213 else if(index1==2 && index2==0 && index3==3 && index4==5){
[5797]214 for(i=0;i<order_horiz;i++){
215 for(j=0;j<order_vert;j++){
[5808]216 coords1[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
217 coords2[i*order_vert+j]=0.0;
218 coords3[i*order_vert+j]= 0.5*(1-seg_horiz_coords[i]);
219 coords4[i*order_vert+j]=seg_vert_coords[j];
220 weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
[5797]221 }
222 }
223 }
224 else{
[12706]225 _error2_("Tria not supported yet (user provided indices " << index1 << " " << index2 << " " << index3 << " " << index4 << ")");
[5797]226 }
227
[5808]228 /*clean-up*/
[12706]229 xDelete<double>(seg_horiz_coords);
230 xDelete<double>(seg_horiz_weights);
231 xDelete<double>(seg_vert_coords);
232 xDelete<double>(seg_vert_weights);
[5797]233}
234/*}}}*/
[12706]235/*FUNCTION GaussPenta::~GaussPenta(){{{*/
[5641]236GaussPenta::~GaussPenta(){
[12706]237 xDelete<double>(weights);
238 xDelete<double>(coords1);
239 xDelete<double>(coords2);
240 xDelete<double>(coords3);
241 xDelete<double>(coords4);
[5641]242}
243/*}}}*/
244
245/*Methods*/
[12706]246/*FUNCTION GaussPenta::Echo{{{*/
[5641]247void GaussPenta::Echo(void){
248
[12706]249 _printLine_("GaussPenta:");
250 _printLine_(" numgauss: " << numgauss);
[5641]251
252 if (weights){
[12706]253 _printString_(" weights = [");
254 for(int i=0;i<numgauss;i++) _printLine_(" " << weights[i]);
255 _printLine_("]");
[5641]256 }
[12706]257 else _printLine_("weights = NULL");
[5641]258 if (coords1){
[12706]259 _printString_(" coords1 = [");
260 for(int i=0;i<numgauss;i++) _printLine_(" " << coords1[i]);
261 _printLine_("]");
[5641]262 }
[12706]263 else _printLine_("coords1 = NULL");
[5641]264 if (coords2){
[12706]265 _printString_(" coords2 = [");
266 for(int i=0;i<numgauss;i++) _printLine_(" " << coords2[i]);
267 _printLine_("]");
[5641]268 }
[12706]269 else _printLine_("coords2 = NULL");
[5641]270 if (coords3){
[12706]271 _printString_(" coords3 = [");
272 for(int i=0;i<numgauss;i++) _printLine_(" " << coords3[i]);
273 _printLine_("]");
[5641]274 }
[12706]275 else _printLine_("coords3 = NULL");
[5641]276 if (coords4){
[12706]277 _printString_(" coords4 = [");
278 for(int i=0;i<numgauss;i++) _printLine_(" " << coords4[i]);
279 _printLine_("]");
[5641]280 }
[12706]281 else _printLine_("coords4 = NULL");
[5641]282
[12706]283 _printLine_(" weight = " << weight);
284 _printLine_(" coord1 = " << coord1);
285 _printLine_(" coord2 = " << coord2);
286 _printLine_(" coord3 = " << coord3);
287 _printLine_(" coord4 = " << coord4);
[5641]288
289}
290/*}}}*/
[12706]291/*FUNCTION GaussPenta::GaussCenter{{{*/
[5641]292void GaussPenta::GaussCenter(void){
293
294 /*update static arrays*/
295 coord1=ONETHIRD;
296 coord2=ONETHIRD;
297 coord3=ONETHIRD;
298 coord4=0.0;
299
300}
301/*}}}*/
[12706]302/*FUNCTION GaussPenta::GaussPoint{{{*/
[5641]303void GaussPenta::GaussPoint(int ig){
304
305 /*Check input in debugging mode*/
[6412]306 _assert_(ig>=0 && ig< numgauss);
[5641]307
308 /*update static arrays*/
309 weight=weights[ig];
310 coord1=coords1[ig];
311 coord2=coords2[ig];
312 coord3=coords3[ig];
313 coord4=coords4[ig];
314
315}
316/*}}}*/
[12706]317/*FUNCTION GaussPenta::GaussVertex{{{*/
[5641]318void GaussPenta::GaussVertex(int iv){
319
320 /*in debugging mode: check that the default constructor has been called*/
[6412]321 _assert_(numgauss==-1);
[5641]322
323 /*update static arrays*/
324 switch(iv){
325 case 0:
326 coord1=1; coord2=0; coord3=0; coord4= -1;
327 break;
328 case 1:
329 coord1=0; coord2=1; coord3=0; coord4= -1;
330 break;
331 case 2:
332 coord1=0; coord2=0; coord3=1; coord4= -1;
333 break;
334 case 3:
335 coord1=1; coord2=0; coord3=0; coord4= +1;
336 break;
337 case 4:
338 coord1=0; coord2=1; coord3=0; coord4= +1;
339 break;
340 case 5:
341 coord1=0; coord2=0; coord3=1; coord4= +1;
342 break;
343 default:
[12706]344 _error2_("vertex index should be in [0 5]");
[5641]345
346 }
347
348}
349/*}}}*/
[12706]350/*FUNCTION GaussPenta::GaussFaceTria{{{*/
[5675]351void GaussPenta::GaussFaceTria(int index1, int index2, int index3, int order){
352
353 /*in debugging mode: check that the default constructor has been called*/
[6412]354 _assert_(numgauss==-1);
[5675]355
356 /*Basal Tria*/
357 if(index1==0 && index2==1 && index3==2){
358 GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
[12706]359 coords4=xNew<double>(numgauss);
[5675]360 for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
361 }
362 else{
[12706]363 _error2_("Tria not supported yet");
[5675]364 }
365
366}
367/*}}}*/
[12706]368/*FUNCTION GaussPenta::begin{{{*/
[5641]369int GaussPenta::begin(void){
370
371 /*Check that this has been initialized*/
[6412]372 _assert_(numgauss>0);
373 _assert_(weights);
374 _assert_(coords1);
375 _assert_(coords2);
376 _assert_(coords3);
377 _assert_(coords4);
[5641]378
379 /*return first gauss index*/
380 return 0;
381}
382/*}}}*/
[12706]383/*FUNCTION GaussPenta::end{{{*/
[5641]384int GaussPenta::end(void){
385
386 /*Check that this has been initialized*/
[6412]387 _assert_(numgauss>0);
388 _assert_(weights);
389 _assert_(coords1);
390 _assert_(coords2);
391 _assert_(coords3);
392 _assert_(coords4);
[5641]393
394 /*return last gauss index +1*/
395 return numgauss;
396}
397/*}}}*/
[12706]398/*FUNCTION GaussPenta::SynchronizeGaussTria{{{*/
[5651]399void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){
400
401 gauss_tria->coord1=this->coord1;
402 gauss_tria->coord2=this->coord2;
403 gauss_tria->coord3=this->coord3;
404 gauss_tria->weight=UNDEF;
405}
406/*}}}*/
Note: See TracBrowser for help on using the repository browser.