source: issm/trunk/src/c/objects/Gauss/GaussTria.cpp@ 5739

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

NumericalFlux does not use beams anymore
File size: 6.3 KB
Line 
1/*!\file GaussTria.c
2 * \brief: implementation of the GaussTria object
3 */
4
5/*Include files: {{{1*/
6#include "./../objects.h"
7/*}}}*/
8
9/*GaussTria constructors and destructors:*/
10/*FUNCTION GaussTria::GaussTria() {{{1*/
11GaussTria::GaussTria(){
12
13 numgauss=-1;
14
15 weights=NULL;
16 coords1=NULL;
17 coords2=NULL;
18 coords3=NULL;
19
20 weight=UNDEF;
21 coord1=UNDEF;
22 coord2=UNDEF;
23 coord3=UNDEF;
24}
25/*}}}*/
26/*FUNCTION GaussTria::GaussTria(int order) {{{1*/
27GaussTria::GaussTria(int order){
28
29 /*Get gauss points*/
30 GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
31
32 /*Initialize static fields as undefinite*/
33 weight=UNDEF;
34 coord1=UNDEF;
35 coord2=UNDEF;
36 coord3=UNDEF;
37
38}
39/*}}}*/
40/*FUNCTION GaussTria::GaussTria(int index1,int index2,int order) {{{1*/
41GaussTria::GaussTria(int index1,int index2,int order){
42
43 /*Intermediaties*/
44 double *seg_coords = NULL;
45 double *seg_weights = NULL;
46 int i,index3;
47
48 /*Get Segment gauss points*/
49 numgauss=order;
50 GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);
51
52 /*Allocate GaussTria fields*/
53 coords1=(double*)xmalloc(numgauss*sizeof(double));
54 coords2=(double*)xmalloc(numgauss*sizeof(double));
55 coords3=(double*)xmalloc(numgauss*sizeof(double));
56 weights=(double*)xmalloc(numgauss*sizeof(double));
57
58 /*Reverse index1 and 2 if necessary*/
59 if (index1>index2){
60 index3=index1; index1=index2; index2=index3;
61 for(i=0;i<numgauss;i++) seg_coords[i]=-seg_coords[i];
62 }
63
64 /*Build Triangle Gauss point*/
65 if (index1==0 && index2==1){
66 for(i=0;i<numgauss;i++) coords1[i]= 0.5*(1-seg_coords[i]);
67 for(i=0;i<numgauss;i++) coords2[i]=1-0.5*(1.-seg_coords[i]);
68 for(i=0;i<numgauss;i++) coords3[i]=0;
69 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
70 }
71 else if (index1==0 && index2==2){
72 for(i=0;i<numgauss;i++) coords1[i]= 0.5*(1-seg_coords[i]);
73 for(i=0;i<numgauss;i++) coords2[i]=0;
74 for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
75 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
76 }
77 else if (index1==1 && index2==2){
78 for(i=0;i<numgauss;i++) coords1[i]=0;
79 for(i=0;i<numgauss;i++) coords2[i]= 0.5*(1-seg_coords[i]);
80 for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
81 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
82 }
83 else
84 ISSMERROR("The 2 indices provided are not supported yet (user provided %i and %i)",index1,index2);
85
86 /*Initialize static fields as undefined*/
87 weight=UNDEF;
88 coord1=UNDEF;
89 coord2=UNDEF;
90 coord3=UNDEF;
91
92 /*clean up*/
93 xfree((void**)&seg_coords);
94 xfree((void**)&seg_weights);
95}
96/*}}}*/
97/*FUNCTION GaussTria::~GaussTria(){{{1*/
98GaussTria::~GaussTria(){
99 xfree((void**)&weights);
100 xfree((void**)&coords1);
101 xfree((void**)&coords2);
102 xfree((void**)&coords3);
103}
104/*}}}*/
105
106/*Methods*/
107/*FUNCTION GaussTria::Echo{{{1*/
108void GaussTria::Echo(void){
109
110 printf("GaussTria:\n");
111 printf(" numgauss: %i\n",numgauss);
112
113 if (weights){
114 printf(" weights = [");
115 for(int i=0;i<numgauss;i++) printf(" %g\n",weights[i]);
116 printf("]\n");
117 }
118 else printf("weights = NULL\n");
119 if (coords1){
120 printf(" coords1 = [");
121 for(int i=0;i<numgauss;i++) printf(" %g\n",coords1[i]);
122 printf("]\n");
123 }
124 else printf("coords1 = NULL\n");
125 if (coords2){
126 printf(" coords2 = [");
127 for(int i=0;i<numgauss;i++) printf(" %g\n",coords2[i]);
128 printf("]\n");
129 }
130 else printf("coords2 = NULL\n");
131 if (coords3){
132 printf(" coords3 = [");
133 for(int i=0;i<numgauss;i++) printf(" %g\n",coords3[i]);
134 printf("]\n");
135 }
136 else printf("coords3 = NULL\n");
137
138 printf(" weight = %g\n",weight);
139 printf(" coord1 = %g\n",coord1);
140 printf(" coord2 = %g\n",coord2);
141 printf(" coord3 = %g\n",coord3);
142
143}
144/*}}}*/
145/*FUNCTION GaussTria::GaussCenter{{{1*/
146void GaussTria::GaussCenter(void){
147
148 /*update static arrays*/
149 coord1=ONETHIRD;
150 coord2=ONETHIRD;
151 coord3=ONETHIRD;
152
153}
154/*}}}*/
155/*FUNCTION GaussTria::GaussEdgeCenter{{{1*/
156void GaussTria::GaussEdgeCenter(int index1,int index2){
157
158 int index3;
159
160 /*Reverse index1 and 2 if necessary*/
161 if (index1>index2){
162 index3=index1; index1=index2; index2=index3;
163 }
164
165 /*update static arrays*/
166 if (index1==0 && index2==1){
167 coord1=0.5;
168 coord2=0.5;
169 coord3=0.0;
170 }
171 else if (index1==0 && index2==2){
172 coord1=0.5;
173 coord2=0.0;
174 coord3=0.5;
175 }
176 else if (index1==1 && index2==2){
177 coord1=0.0;
178 coord2=0.5;
179 coord3=0.5;
180 }
181 else
182 ISSMERROR("The 2 indices provided are not supported yet (user provided %i and %i)",index1,index2);
183
184}
185/*}}}*/
186/*FUNCTION GaussTria::GaussPoint{{{1*/
187void GaussTria::GaussPoint(int ig){
188
189 /*Check input in debugging mode*/
190 ISSMASSERT(ig>=0 && ig< numgauss);
191
192 /*update static arrays*/
193 weight=weights[ig];
194 coord1=coords1[ig];
195 coord2=coords2[ig];
196 coord3=coords3[ig];
197
198}
199/*}}}*/
200/*FUNCTION GaussTria::GaussFromCoords{{{1*/
201void GaussTria::GaussFromCoords(double x,double y,double* xyz_list){
202
203 /*Intermediaries*/
204 double area = 0;
205 double x1,y1,x2,y2,x3,y3;
206
207 /*in debugging mode: check that the default constructor has been called*/
208 ISSMASSERT(numgauss==-1);
209
210 x1=*(xyz_list+3*0+0); y1=*(xyz_list+3*0+1);
211 x2=*(xyz_list+3*1+0); y2=*(xyz_list+3*1+1);
212 x3=*(xyz_list+3*2+0); y3=*(xyz_list+3*2+1);
213
214 area=(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
215
216 /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/
217 coord1=((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area;
218
219 /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/
220 coord2=((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;
221
222 /*Get third area coordinate 1-area1-area2: */
223 coord3=1-coord1-coord2;
224
225}
226/*}}}*/
227/*FUNCTION GaussTria::GaussVertex{{{1*/
228void GaussTria::GaussVertex(int iv){
229
230 /*in debugging mode: check that the default constructor has been called*/
231 ISSMASSERT(numgauss==-1);
232
233 /*update static arrays*/
234 switch(iv){
235 case 0:
236 coord1=1; coord2=0; coord3=0;
237 break;
238 case 1:
239 coord1=0; coord2=1; coord3=0;
240 break;
241 case 2:
242 coord1=0; coord2=0; coord3=1;
243 break;
244 default:
245 ISSMERROR("vertex index should be in [0 2]");
246
247 }
248
249}
250/*}}}*/
251/*FUNCTION GaussTria::begin{{{1*/
252int GaussTria::begin(void){
253
254 /*Check that this has been initialized*/
255 ISSMASSERT(numgauss>0);
256 ISSMASSERT(weights);
257 ISSMASSERT(coords1);
258 ISSMASSERT(coords2);
259 ISSMASSERT(coords3);
260
261 /*return first gauss index*/
262 return 0;
263}
264/*}}}*/
265/*FUNCTION GaussTria::end{{{1*/
266int GaussTria::end(void){
267
268 /*Check that this has been initialized*/
269 ISSMASSERT(numgauss>0);
270 ISSMASSERT(weights);
271 ISSMASSERT(coords1);
272 ISSMASSERT(coords2);
273 ISSMASSERT(coords3);
274
275 /*return last gauss index +1*/
276 return numgauss;
277}
278/*}}}*/
Note: See TracBrowser for help on using the repository browser.