xref: /phasta/phSolver/common/symquadw.c (revision 97a07b0ae3d9fb83b0ec5c51c6b02ee38dcf6667)
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 
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 
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