1 #include <stdlib.h>
2
3 #include <FCMangle.h>
4 #define symquad FortranCInterface_GLOBAL_(symquad, SYMQUAD)
5
6 typedef double DARR3[3];
7
8 int quadIntPnt(int,DARR3**,double **);
9
symquad(int * n1,double pt[][4],double wt[],int * err)10 void symquad(int *n1, double pt[][4], double wt[], int *err)
11 {
12 double *lwt;
13 DARR3 *lpt;
14 int i,j;
15 *err = quadIntPnt(*n1, &lpt, &lwt);
16 for(i=0; i < *n1; i++) {
17 wt[i] = lwt[i];
18 for(j=0; j < 3; j++)
19 pt[i][j]=lpt[i][j];
20 }
21 }
22
23 /*$Id$*/
24 #include <stdio.h>
25
26 #define QpNon -1.000000000000000
27 #define Qp21 -0.577350269189626
28 #define Qp22 0.577350269189626
29 #define Qw2 1.000000000000000
30
31 #define Qp31 -0.774596669241483
32 #define Qp32 0.000000000000000
33 #define Qp33 0.774596669241483
34 #define Qw31 0.555555555555556
35 #define Qw32 0.888888888888889
36 #define Qw33 0.555555555555556
37 #define Qw311 0.308641975308642
38 #define Qw321 0.493827160493831
39 #define Qw322 0.790123456790124
40 #define Qw312 0.493827160493831
41 #define Qw313 0.308641975308642
42 #define Qw323 0.493827160493831
43 #define Qw331 0.308641975308642
44 #define Qw332 0.493827160493831
45 #define Qw333 0.308641975308642
46
47
48 static double rstw4[][3] = {
49 {Qp21, Qp21, QpNon},
50 {Qp22, Qp21, QpNon},
51 {Qp21, Qp22, QpNon},
52 {Qp22, Qp22, QpNon}
53 };
54
55 static double twt4[] = {Qw2,Qw2,Qw2,Qw2};
56
57 static double rstw9[][3] = {
58 { Qp31, Qp31, QpNon},
59 { Qp32, Qp31, QpNon},
60 { Qp33, Qp31, QpNon},
61 { Qp31, Qp32, QpNon},
62 { Qp32, Qp32, QpNon},
63 { Qp33, Qp32, QpNon},
64 { Qp31, Qp33, QpNon},
65 { Qp32, Qp33, QpNon},
66 { Qp33, Qp33, QpNon}
67 };
68
69 static double twt9[] = { Qw311, Qw321, Qw331, Qw312, Qw322, Qw332, Qw313,
70 Qw323, Qw333 };
71
72
73 #ifdef __cplusplus
74 extern "C" {
75 #endif
76
quadIntPnt(int nint,DARR3 ** bcord,double ** wt)77 int quadIntPnt(int nint, DARR3 **bcord, double **wt)
78 {
79 int retval = 1;
80 int i,j,l;
81 DARR3* rstw;
82 double* twt;
83
84 /* Rule 4 & 5*/
85 /* rule 5 has not been tested */
86 double bp4[]={-0.8611363115940526,-0.3399810435848563,0.339981043584856,0.8611363115940526};
87 double bp5[]={-0.9061798459386640,-0.5384693101056831,0,0.5384693101056831,0.9061798459386640};
88 double bw4[]={0.3478548451374539,0.6521451548625461,0.6521451548625461,0.3478548451374539};
89 double bw5[]={0.2369268850561891,0.4786286704993665,0.5688888888888889,0.4786286704993665,0.2369268850561891};
90
91 if( nint == 4 ){*bcord = rstw4 ; *wt = twt4; }
92 else if (nint == 9){*bcord = rstw9 ; *wt = twt9; }
93 else if (nint == 16){
94
95 rstw = (DARR3 *)malloc(16*sizeof(DARR3));
96 twt = (double *)malloc(16*sizeof(double));
97
98 i=j=0;
99 for(l=0;l<16;l++){
100 rstw[l][0]=bp4[i];
101 rstw[l][1]=bp4[j];
102 rstw[l][2]=QpNon;
103 twt[l]= bw4[i]*bw4[j];
104 i++;
105 if( i == 4) { i=0; j++;}
106 }
107 *bcord = rstw; *wt = twt;
108
109 }else if(nint == 25){
110
111 rstw = (DARR3 *)malloc(25*sizeof(DARR3));
112 twt = (double *)malloc(25*sizeof(double));
113
114 i=j=0;
115 for(l=0;l<25;l++){
116 rstw[l][0]=bp5[i];
117 rstw[l][1]=bp5[j];
118 rstw[l][2]=QpNon;
119 twt[l]= bw5[i]*bw5[j];
120 i++;
121 if( i == 5) { i=0; j++;}
122 }
123 *bcord = rstw; *wt = twt;
124
125 } else {
126
127 fprintf(stderr,"\n%d integration points unsupported in symquad.c; give {4,9,16,25}\n",nint);
128 retval = 0;
129 }
130 return retval ;
131 }
132
133 #ifdef __cplusplus
134 }
135 #endif
136