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*/
|
---|
11 | GaussTria::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*/
|
---|
27 | GaussTria::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*/
|
---|
41 | GaussTria::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*/
|
---|
98 | GaussTria::~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*/
|
---|
108 | void 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*/
|
---|
146 | void 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*/
|
---|
156 | void 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*/
|
---|
187 | void 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*/
|
---|
201 | void 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*/
|
---|
228 | void 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*/
|
---|
252 | int 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*/
|
---|
266 | int 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 | /*}}}*/
|
---|