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 grid1 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::GaussPoint{{{1*/
|
---|
156 | void GaussTria::GaussPoint(int ig){
|
---|
157 |
|
---|
158 | /*Check input in debugging mode*/
|
---|
159 | ISSMASSERT(ig>=0 && ig< numgauss);
|
---|
160 |
|
---|
161 | /*update static arrays*/
|
---|
162 | weight=weights[ig];
|
---|
163 | coord1=coords1[ig];
|
---|
164 | coord2=coords2[ig];
|
---|
165 | coord3=coords3[ig];
|
---|
166 |
|
---|
167 | }
|
---|
168 | /*}}}*/
|
---|
169 | /*FUNCTION GaussTria::GaussVertex{{{1*/
|
---|
170 | void GaussTria::GaussVertex(int iv){
|
---|
171 |
|
---|
172 | /*in debugging mode: check that the default constructor has been called*/
|
---|
173 | ISSMASSERT(numgauss==-1);
|
---|
174 |
|
---|
175 | /*update static arrays*/
|
---|
176 | switch(iv){
|
---|
177 | case 0:
|
---|
178 | coord1=1; coord2=0; coord3=0;
|
---|
179 | break;
|
---|
180 | case 1:
|
---|
181 | coord1=0; coord2=1; coord3=0;
|
---|
182 | break;
|
---|
183 | case 2:
|
---|
184 | coord1=0; coord2=0; coord3=1;
|
---|
185 | break;
|
---|
186 | default:
|
---|
187 | ISSMERROR("vertex index should be in [0 2]");
|
---|
188 |
|
---|
189 | }
|
---|
190 |
|
---|
191 | }
|
---|
192 | /*}}}*/
|
---|
193 | /*FUNCTION GaussTria::begin{{{1*/
|
---|
194 | int GaussTria::begin(void){
|
---|
195 |
|
---|
196 | /*Check that this has been initialized*/
|
---|
197 | ISSMASSERT(numgauss>0);
|
---|
198 | ISSMASSERT(weights);
|
---|
199 | ISSMASSERT(coords1);
|
---|
200 | ISSMASSERT(coords2);
|
---|
201 | ISSMASSERT(coords3);
|
---|
202 |
|
---|
203 | /*return first gauss index*/
|
---|
204 | return 0;
|
---|
205 | }
|
---|
206 | /*}}}*/
|
---|
207 | /*FUNCTION GaussTria::end{{{1*/
|
---|
208 | int GaussTria::end(void){
|
---|
209 |
|
---|
210 | /*Check that this has been initialized*/
|
---|
211 | ISSMASSERT(numgauss>0);
|
---|
212 | ISSMASSERT(weights);
|
---|
213 | ISSMASSERT(coords1);
|
---|
214 | ISSMASSERT(coords2);
|
---|
215 | ISSMASSERT(coords3);
|
---|
216 |
|
---|
217 | /*return last gauss index +1*/
|
---|
218 | return numgauss;
|
---|
219 | }
|
---|
220 | /*}}}*/
|
---|