1 #include <FCMangle.h>
2 #define symquadw FortranCInterface_GLOBAL_(symquadw,SYMQUADW)
3
4 typedef double DARR3[3];
5
6 int quadwIntPnt(int,DARR3**,double **);
7
symquadw(int * n1,double pt[][4],double wt[],int * err)8 void symquadw(int *n1, double pt[][4], double wt[], int *err)
9 {
10 double *lwt = 0;
11 DARR3 *lpt = 0;
12 int i,j;
13 *err = quadwIntPnt(*n1, &lpt, &lwt);
14 for(i=0; i < *n1; i++) {
15 wt[i] = lwt[i];
16 for(j=0; j < 3; j++)
17 pt[i][j]=lpt[i][j];
18 }
19 }
20
21 /*$Id$*/
22 #include <stdio.h>
23
24 /* Rule 2 constants */
25
26 #define Qp21 -0.577350269189626
27 #define Qp22 0.577350269189626
28 #define Pp21 0.211324865405187
29 #define Pp22 0.788675134594813
30 #define Qw2 1.000000000000000
31
32 /* Rule 3 constants */
33
34 #define Qp31 -0.774596669241483
35 #define Qp32 0.000000000000000
36 #define Qp33 0.774596669241483
37 #define Qw311 0.308641975308642 /* Qw31 * Qw31 note: Qw31=Qw33*/
38 #define Qw321 0.493827160493831 /* Qw32 * Qw31 */
39 #define Qw322 0.790123456790124 /* Qw32 * Qw32 */
40 #define Pp31 0.112701665379259
41 #define Pp32 0.500000000000000
42 #define Pp33 0.887298334620741
43
44 /* Rule 4 constants */
45
46 #define Qp41 -0.861136311594053
47 #define Qp42 -0.339981043584856
48 #define Qp43 0.339981043584856
49 #define Qp44 0.861136311594053
50 #define Qw41 0.347854845137454
51 #define Qw42 0.652145154862544
52 #define Qw43 0.652145154862544
53 #define Qw44 0.347854845137544
54 #define Pp41 0.873821971016996
55 #define Pp42 0.063089014491502
56 #define Pp43 0.501426509658179
57 #define Pp44 0.249286745170910
58 #define Qw4141 0.121002993285602 /* Qw41 * Qw41 */
59 #define Qw4142 0.226851851851851 /* Qw41 * Qw42 */
60 #define Qw4143 0.226851851851851 /* Qw41 * Qw43 */
61 #define Qw4144 0.121002993285602 /* Qw41 * Qw44 */
62 #define Qw4241 0.226851851851851 /* etc.. */
63 #define Qw4242 0.425293303010692
64 #define Qw4243 0.425293303010692
65 #define Qw4244 0.226851851851851
66 #define Qw4341 0.226851851851851
67 #define Qw4342 0.425293303010692
68 #define Qw4343 0.425293303010692
69 #define Qw4344 0.226851851851851
70 #define Qw4441 0.121002993285602
71 #define Qw4442 0.226851851851851
72 #define Qw4443 0.226851851851851
73 #define Qw4444 0.121002993285602
74
75
76 /* typedef double DARR3[3] ; */
77
78 /* Rule 2 */
79
80 static double rstw4[][3] = {
81 {Pp21, 0.0 , Qp21},
82 {Pp22, 0.0 , Qp21},
83 {Pp21, 0.0 , Qp22},
84 {Pp22, 0.0 , Qp22}
85 };
86
87 static double twt4[] = {Qw2,Qw2,Qw2,Qw2};
88
89 /* Rule 3 */
90
91 static double rstw9[][3] = {
92 {Pp31,0.0,Qp31},
93 {Pp32,0.0,Qp31},
94 {Pp33,0.0,Qp31},
95 {Pp31,0.0,Qp32},
96 {Pp32,0.0,Qp32},
97 {Pp33,0.0,Qp32},
98 {Pp31,0.0,Qp33},
99 {Pp32,0.0,Qp33},
100 {Pp33,0.0,Qp33}
101 };
102
103 static double twt9[] =
104 {Qw311, Qw321, Qw311, Qw321, Qw322, Qw321, Qw311, Qw321, Qw311};
105
106 /* Rule 4 */
107
108 static double rstw16[][3] = {
109 {Pp41,0.0,Qp41},
110 {Pp42,0.0,Qp41},
111 {Pp43,0.0,Qp41},
112 {Pp44,0.0,Qp41},
113 {Pp41,0.0,Qp42},
114 {Pp42,0.0,Qp42},
115 {Pp43,0.0,Qp42},
116 {Pp44,0.0,Qp42},
117 {Pp41,0.0,Qp43},
118 {Pp42,0.0,Qp43},
119 {Pp43,0.0,Qp43},
120 {Pp44,0.0,Qp43},
121 {Pp41,0.0,Qp44},
122 {Pp42,0.0,Qp44},
123 {Pp43,0.0,Qp44},
124 {Pp44,0.0,Qp44}
125 };
126
127 static double twt16[] =
128 {Qw4141, Qw4241, Qw4341, Qw4441, Qw4142, Qw4242, Qw4342, Qw4442,
129 Qw4143, Qw4243, Qw4343, Qw4443, Qw4144, Qw4244, Qw4344, Qw4444};
130
131
132 #ifdef __cplusplus
133 extern "C" {
134 #endif
135
quadwIntPnt(int nint,DARR3 ** bcord,double ** wt)136 int quadwIntPnt(int nint, DARR3 **bcord, double **wt)
137 {
138 int retval = 1;
139
140 if( nint == 4 ){*bcord = rstw4 ; *wt = twt4; }
141 else if( nint == 9 ){*bcord = rstw9 ; *wt = twt9; }
142 else if( nint == 16){*bcord = rstw16 ; *wt = twt16;}
143 else
144 {
145 fprintf(stderr,"\n%d integration points unsupported in symquadw.c; give {4,9,16}\n",nint);
146 retval = 0;
147 }
148 return retval ;
149 }
150
151 #ifdef __cplusplus
152 }
153 #endif
154