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

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

Deleted some functions and added more GaussPenta
File size: 6.4 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{{{1*/
83GaussPenta::GaussPenta(int index1, int index2, int index3, int order){
84
85 /*Basal Tria*/
86 if(index1==0 && index2==1 && index3==2){
87
88 /*Get GaussTria*/
89 GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
90
91 /*compute z coordinate*/
92 coords4=(double*)xmalloc(numgauss*sizeof(double));
93 for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
94 }
95 else{
96 ISSMERROR("Tria not supported yet");
97 }
98
99}
100/*}}}*/
101/*FUNCTION GaussPenta::~GaussPenta(){{{1*/
102GaussPenta::~GaussPenta(){
103 xfree((void**)&weights);
104 xfree((void**)&coords1);
105 xfree((void**)&coords2);
106 xfree((void**)&coords3);
107 xfree((void**)&coords4);
108}
109/*}}}*/
110
111/*Methods*/
112/*FUNCTION GaussPenta::Echo{{{1*/
113void GaussPenta::Echo(void){
114
115 printf("GaussPenta:\n");
116 printf(" numgauss: %i\n",numgauss);
117
118 if (weights){
119 printf(" weights = [");
120 for(int i=0;i<numgauss;i++) printf(" %g\n",weights[i]);
121 printf("]\n");
122 }
123 else printf("weights = NULL\n");
124 if (coords1){
125 printf(" coords1 = [");
126 for(int i=0;i<numgauss;i++) printf(" %g\n",coords1[i]);
127 printf("]\n");
128 }
129 else printf("coords1 = NULL\n");
130 if (coords2){
131 printf(" coords2 = [");
132 for(int i=0;i<numgauss;i++) printf(" %g\n",coords2[i]);
133 printf("]\n");
134 }
135 else printf("coords2 = NULL\n");
136 if (coords3){
137 printf(" coords3 = [");
138 for(int i=0;i<numgauss;i++) printf(" %g\n",coords3[i]);
139 printf("]\n");
140 }
141 else printf("coords3 = NULL\n");
142 if (coords4){
143 printf(" coords4 = [");
144 for(int i=0;i<numgauss;i++) printf(" %g\n",coords4[i]);
145 printf("]\n");
146 }
147 else printf("coords4 = NULL\n");
148
149 printf(" weight = %g\n",weight);
150 printf(" coord1 = %g\n",coord1);
151 printf(" coord2 = %g\n",coord2);
152 printf(" coord3 = %g\n",coord3);
153 printf(" coord4 = %g\n",coord4);
154
155}
156/*}}}*/
157/*FUNCTION GaussPenta::GaussCenter{{{1*/
158void GaussPenta::GaussCenter(void){
159
160 /*update static arrays*/
161 coord1=ONETHIRD;
162 coord2=ONETHIRD;
163 coord3=ONETHIRD;
164 coord4=0.0;
165
166}
167/*}}}*/
168/*FUNCTION GaussPenta::GaussPoint{{{1*/
169void GaussPenta::GaussPoint(int ig){
170
171 /*Check input in debugging mode*/
172 ISSMASSERT(ig>=0 && ig< numgauss);
173
174 /*update static arrays*/
175 weight=weights[ig];
176 coord1=coords1[ig];
177 coord2=coords2[ig];
178 coord3=coords3[ig];
179 coord4=coords4[ig];
180
181}
182/*}}}*/
183/*FUNCTION GaussPenta::GaussVertex{{{1*/
184void GaussPenta::GaussVertex(int iv){
185
186 /*in debugging mode: check that the default constructor has been called*/
187 ISSMASSERT(numgauss==-1);
188
189 /*update static arrays*/
190 switch(iv){
191 case 0:
192 coord1=1; coord2=0; coord3=0; coord4= -1;
193 break;
194 case 1:
195 coord1=0; coord2=1; coord3=0; coord4= -1;
196 break;
197 case 2:
198 coord1=0; coord2=0; coord3=1; coord4= -1;
199 break;
200 case 3:
201 coord1=1; coord2=0; coord3=0; coord4= +1;
202 break;
203 case 4:
204 coord1=0; coord2=1; coord3=0; coord4= +1;
205 break;
206 case 5:
207 coord1=0; coord2=0; coord3=1; coord4= +1;
208 break;
209 default:
210 ISSMERROR("vertex index should be in [0 5]");
211
212 }
213
214}
215/*}}}*/
216/*FUNCTION GaussPenta::GaussFaceTria{{{1*/
217void GaussPenta::GaussFaceTria(int index1, int index2, int index3, int order){
218
219 /*in debugging mode: check that the default constructor has been called*/
220 ISSMASSERT(numgauss==-1);
221
222 /*Basal Tria*/
223 if(index1==0 && index2==1 && index3==2){
224 GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
225 coords4=(double*)xmalloc(numgauss*sizeof(double));
226 for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
227 }
228 else{
229 ISSMERROR("Tria not supported yet");
230 }
231
232}
233/*}}}*/
234/*FUNCTION GaussPenta::begin{{{1*/
235int GaussPenta::begin(void){
236
237 /*Check that this has been initialized*/
238 ISSMASSERT(numgauss>0);
239 ISSMASSERT(weights);
240 ISSMASSERT(coords1);
241 ISSMASSERT(coords2);
242 ISSMASSERT(coords3);
243 ISSMASSERT(coords4);
244
245 /*return first gauss index*/
246 return 0;
247}
248/*}}}*/
249/*FUNCTION GaussPenta::end{{{1*/
250int GaussPenta::end(void){
251
252 /*Check that this has been initialized*/
253 ISSMASSERT(numgauss>0);
254 ISSMASSERT(weights);
255 ISSMASSERT(coords1);
256 ISSMASSERT(coords2);
257 ISSMASSERT(coords3);
258 ISSMASSERT(coords4);
259
260 /*return last gauss index +1*/
261 return numgauss;
262}
263/*}}}*/
264/*FUNCTION GaussPenta::SynchronizeGaussTria{{{1*/
265void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){
266
267 gauss_tria->coord1=this->coord1;
268 gauss_tria->coord2=this->coord2;
269 gauss_tria->coord3=this->coord3;
270 gauss_tria->weight=UNDEF;
271}
272/*}}}*/
Note: See TracBrowser for help on using the repository browser.