1 #include <FCMangle.h>
2 #define symwdg FortranCInterface_GLOBAL_(symwdg,SYMWDG)
3
4 typedef double DARR4[4];
5
6 int wdgIntPnt(int,DARR4**,double **);
7
symwdg(int * n1,double pt[][4],double wt[],int * err)8 void symwdg(int *n1, double pt[][4], double wt[], int *err)
9 {
10 double *lwt;
11 DARR4 *lpt;
12 int i,j;
13 *err = wdgIntPnt(*n1, &lpt, &lwt);
14 for(i=0; i < *n1; i++) {
15 wt[i] = lwt[i];
16 for(j=0; j <4; j++)
17 pt[i][j]=lpt[i][j];
18 }
19 }
20
21 /*$Id$*/
22 #include <stdio.h>
23
24 /* constants for 3D - 2pt rule */
25
26 #define Qp21 -0.577350269189626
27 #define Qp22 0.577350269189626
28 #define Pp23 0.166666666666667
29 #define Pp24 0.666666666666667
30 #define Ww32 0.666666666666667 /* Pw22 * Qw2 */
31
32 /* constants for 3D - rule 3 (18 point) */
33
34 #define Qp31 -0.774596669241483
35 #define Qp32 0.000000000000000
36 #define Qp33 0.774596669241483
37 #define Qw31 0.555555555555556
38 #define Qw32 0.888888888888889
39 #define Qw33 0.555555555555556
40
41 #define Pp34 0.816847572980459
42 #define Pp35 0.091576213509771
43 #define Pp36 0.108103018168070
44 #define Pp37 0.445948490915965
45 #define Pw33 0.219903487310644
46 #define Pw34 0.446763179356022
47
48 #define Ww3331 0.122168604061469 /* Pw33 * Qw31 */
49 #define Ww3332 0.195469766498350 /* Pw33 * Qw32 */
50 #define Ww3333 0.122168604061469 /* Pw33 * Qw33 */
51 #define Ww3431 0.248201766308901 /* Pw34 * Qw31 */
52 #define Ww3432 0.397122826094242 /* Pw34 * Qw32 */
53 #define Ww3433 0.248201766308901 /* Pw34 * Qw33 */
54
55 /* constants for 3D rule 4 (48 point) */
56
57 #define Qp41 -0.861136311594053
58 #define Qp42 -0.339981043584856
59 #define Qp43 0.339981043584856
60 #define Qp44 0.861136311594053
61 #define Qw41 0.347854845137454
62 #define Qw42 0.652145154862544
63 #define Qw43 0.652145154862544
64 #define Qw44 0.347854845137544
65
66 #define Pp41 0.873821971016996
67 #define Pp42 0.063089014491502
68 #define Pp43 0.501426509658179
69 #define Pp44 0.249286745170910
70 #define Pp45 0.636502499121399
71 #define Pp46 0.310352451033785
72 #define Pp47 0.053145049844816
73 #define Pw41 0.101689812740414 /* S in symtri.c times 2 */
74 #define Pw42 0.233572551452758 /* T in symtri.c times 2 */
75 #define Pw43 0.165702151236748 /* U in symtri.c times 2 */
76
77 /* the weights should go */
78 /*(Pw41,Pw41,Pw41,Pw42,Pw42,Pw42,Pw43,Pw43,Pw43,Pw43,Pw43,Pw43)*Qw41 */
79 /* >> *Qw42,Qw43,Qw44 */
80
81 #define Ww4141 0.0353732940628734 /* Pw41 * Qw41 */
82 #define Ww4142 0.0663165186775404 /* Pw41 * Qw42 */
83 #define Ww4143 0.0663165186775404 /* Pw41 * Qw43 */
84 #define Ww4144 0.0353732940628734 /* Pw41 * Qw44 */
85 #define Ww4241 0.0812493437139591 /* Pw42 * Qw41 */
86 #define Ww4242 0.152323207738798 /* Pw42 * Qw42 */
87 #define Ww4243 0.152323207738798 /* Pw42 * Qw43 */
88 #define Ww4244 0.0812493437139591 /* Pw42 * Qw44 */
89 #define Ww4341 0.057640296157402 /* Pw43 * Qw41 */
90 #define Ww4342 0.108061855079346 /* Pw43 * Qw42 */
91 #define Ww4343 0.108061855079346 /* Pw43 * Qw43 */
92 #define Ww4344 0.057640296157402 /* Pw43 * Qw44 */
93
94 /* typedef double DARR4[4] ; */
95
96 /* Rule 1*/
97
98 static double rstw1[][4] = {
99 { 0.333333333333, 0.333333333333, 0.0, 0.0 }
100 };
101
102 static double twt1[] = { 4.0000000 };
103
104
105
106 /* Rule 2 */
107
108 static double rstw6[][4] = {
109 {Pp23,Pp23,Qp21,0.0},
110 {Pp24,Pp23,Qp21,0.0},
111 {Pp23,Pp24,Qp21,0.0},
112 {Pp23,Pp23,Qp22,0.0},
113 {Pp24,Pp23,Qp22,0.0},
114 {Pp23,Pp24,Qp22,0.0}
115 };
116
117 static double twt6[] = {Ww32,Ww32,Ww32,Ww32,Ww32,Ww32};
118
119 /* Rule 3 */
120
121 static double rstw18[][4] = {
122 {Pp34,Pp35,Qp31,0.0},
123 {Pp35,Pp34,Qp31,0.0},
124 {Pp35,Pp35,Qp31,0.0},
125 {Pp36,Pp37,Qp31,0.0},
126 {Pp37,Pp36,Qp31,0.0},
127 {Pp37,Pp37,Qp31,0.0},
128 {Pp34,Pp35,Qp32,0.0},
129 {Pp35,Pp34,Qp32,0.0},
130 {Pp35,Pp35,Qp32,0.0},
131 {Pp36,Pp37,Qp32,0.0},
132 {Pp37,Pp36,Qp32,0.0},
133 {Pp37,Pp37,Qp32,0.0},
134 {Pp34,Pp35,Qp33,0.0},
135 {Pp35,Pp34,Qp33,0.0},
136 {Pp35,Pp35,Qp33,0.0},
137 {Pp36,Pp37,Qp33,0.0},
138 {Pp37,Pp36,Qp33,0.0},
139 {Pp37,Pp37,Qp33,0.0}
140 };
141
142 static double twt18[] =
143 { Ww3331, Ww3331, Ww3331, Ww3431, Ww3431, Ww3431, Ww3332, Ww3332, Ww3332,
144 Ww3432, Ww3432, Ww3432, Ww3333, Ww3333, Ww3333, Ww3433, Ww3433, Ww3433};
145
146 /* Rule 4 */
147
148 static double rstw48[][4] = {
149 {Pp41,Pp42,Qp41,0.0},
150 {Pp42,Pp41,Qp41,0.0},
151 {Pp42,Pp42,Qp41,0.0},
152 {Pp43,Pp44,Qp41,0.0},
153 {Pp44,Pp43,Qp41,0.0},
154 {Pp44,Pp44,Qp41,0.0},
155 {Pp45,Pp46,Qp41,0.0},
156 {Pp46,Pp45,Qp41,0.0},
157 {Pp45,Pp47,Qp41,0.0},
158 {Pp46,Pp47,Qp41,0.0},
159 {Pp47,Pp46,Qp41,0.0},
160 {Pp47,Pp45,Qp41,0.0},
161 {Pp41,Pp42,Qp42,0.0},
162 {Pp42,Pp41,Qp42,0.0},
163 {Pp42,Pp42,Qp42,0.0},
164 {Pp43,Pp44,Qp42,0.0},
165 {Pp44,Pp43,Qp42,0.0},
166 {Pp44,Pp44,Qp42,0.0},
167 {Pp45,Pp46,Qp42,0.0},
168 {Pp46,Pp45,Qp42,0.0},
169 {Pp45,Pp47,Qp42,0.0},
170 {Pp46,Pp47,Qp42,0.0},
171 {Pp47,Pp46,Qp42,0.0},
172 {Pp47,Pp45,Qp42,0.0},
173 {Pp41,Pp42,Qp43,0.0},
174 {Pp42,Pp41,Qp43,0.0},
175 {Pp42,Pp42,Qp43,0.0},
176 {Pp43,Pp44,Qp43,0.0},
177 {Pp44,Pp43,Qp43,0.0},
178 {Pp44,Pp44,Qp43,0.0},
179 {Pp45,Pp46,Qp43,0.0},
180 {Pp46,Pp45,Qp43,0.0},
181 {Pp45,Pp47,Qp43,0.0},
182 {Pp46,Pp47,Qp43,0.0},
183 {Pp47,Pp46,Qp43,0.0},
184 {Pp47,Pp45,Qp43,0.0},
185 {Pp41,Pp42,Qp44,0.0},
186 {Pp42,Pp41,Qp44,0.0},
187 {Pp42,Pp42,Qp44,0.0},
188 {Pp43,Pp44,Qp44,0.0},
189 {Pp44,Pp43,Qp44,0.0},
190 {Pp44,Pp44,Qp44,0.0},
191 {Pp45,Pp46,Qp44,0.0},
192 {Pp46,Pp45,Qp44,0.0},
193 {Pp45,Pp47,Qp44,0.0},
194 {Pp46,Pp47,Qp44,0.0},
195 {Pp47,Pp46,Qp44,0.0},
196 {Pp47,Pp45,Qp44,0.0}
197 };
198
199 static double twt48[] =
200 { Ww4141, Ww4141, Ww4141, Ww4241, Ww4241, Ww4241, Ww4341, Ww4341, Ww4341,
201 Ww4341, Ww4341, Ww4341, Ww4142, Ww4142, Ww4142, Ww4242, Ww4242, Ww4242,
202 Ww4342, Ww4342, Ww4342, Ww4342, Ww4342, Ww4342, Ww4143, Ww4143, Ww4143,
203 Ww4243, Ww4243, Ww4243, Ww4343, Ww4343, Ww4343, Ww4343, Ww4343, Ww4343,
204 Ww4144, Ww4144, Ww4144, Ww4244, Ww4244, Ww4244, Ww4344, Ww4344, Ww4344,
205 Ww4344, Ww4344, Ww4344};
206
207 #ifdef __cplusplus
208 extern "C" {
209 #endif
210
wdgIntPnt(int nint,DARR4 ** bcord,double ** wt)211 int wdgIntPnt(int nint, DARR4 **bcord, double **wt)
212 {
213 int retval = 1;
214
215 if( nint == 6 ){*bcord = rstw6 ; *wt = twt6; }
216 else if( nint == 1 ){*bcord = rstw1 ; *wt = twt1; }
217 else if( nint == 18 ){*bcord = rstw18 ; *wt = twt18; }
218 else if( nint == 48 ){*bcord = rstw48 ; *wt = twt48; }
219 else
220 {
221 fprintf(stderr,"\n%d integration points unsupported symwdg.c; give {6,18,48}\n",nint);
222 retval = 0;
223 }
224 return retval ;
225 }
226
227 #ifdef __cplusplus
228 }
229 #endif
230