1 #include <FCMangle.h>
2 #define symtripyr FortranCInterface_GLOBAL_(symtripyr, SYMTRIPYR)
3
4 typedef double DARR3[3];
5
6 int triIntPntPyr(int, DARR3**,double**);
7
symtripyr(int * n1,double pt[][4],double wt[],int * err)8 void symtripyr(int *n1, double pt[][4], double wt[], int *err)
9 {
10 double *lwt = 0;
11 DARR3 *lpt = 0;
12 int i,j;
13 *err = triIntPntPyr(*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 #include <stdio.h>
22 /* 1pt constants */
23 #define zero 0.0000000000000000
24 #define ppt500 0.5000000000000000
25 #define p2pt00 2.0000000000000000
26
27 /* 4pt constants */
28 #define ppt455 0.4553418012614798
29 #define ppt788 0.7886751345948130
30 #define ppt577 0.5773502691896260
31 #define ppt122 0.1220084679281462
32 #define ppt211 0.2113248654051870
33
34 /* Rule 3 */
35 #define mpt68 -0.6872983346207413
36 #define pt68 0.6872983346207413
37 #define mpt77 -0.7745966692414831
38 #define pt77 0.7745966692414831
39 #define mpt38 -0.3872983346207416
40 #define pt38 0.3872983346207416
41 #define mpt08 -0.08729833462074179
42 #define pt08 0.08729833462074179
43 #define mpt5 -0.5000000000000000
44
45 #define mpt11 -0.1127016653792584
46 #define mpt88 -0.8872983346207415
47
48 #define Qw311 0.308641975308642
49 #define Qw321 0.493827160493831
50 #define Qw322 0.790123456790124
51 #define Qw312 0.493827160493831
52 #define Qw313 0.308641975308642
53 #define Qw323 0.493827160493831
54 #define Qw331 0.308641975308642
55 #define Qw332 0.493827160493831
56 #define Qw333 0.308641975308642
57
58 /* #define J 0.108103018168070 */
59 /* #define K 0.445948490915965 */
60
61 /* #define L 0.225000000000000 */
62 /* #define M 0.125939180544827 */
63 /* #define N 0.132394152788506 */
64
65 /* #define O 0.797426985353087 */
66 /* #define P 0.101286507323456 */
67 /* #define Q 0.470142064105115 */
68 /* #define R 0.059715871789770 */
69
70 /* #define S 0.050844906370207 */
71 /* #define T 0.116786275726379 */
72 /* #define U 0.082851075618374 */
73
74 /* #define V 0.873821971016996 */
75 /* #define W 0.063089014491502 */
76 /* #define X 0.501426509658179 */
77 /* #define Y 0.249286745170910 */
78 /* #define Z 0.636502499121399 */
79 /* #define AA 0.310352451033785 */
80 /* #define BB 0.053145049844816 */
81
82 /* typedef double DARR3[3] ; */
83
84 static double rst1[][3] = {{zero,-1,zero}};
85 static double wt1[] = {2.2360679774997897};
86 /* only rule 2 is implemented correctly at this time */
87 static double rst4[][3] ={
88 {-.455341801261479547, -.788675134594812883 ,-.577350269189625763},
89 { .455341801261479547, -.788675134594812883 ,-.577350269189625763},
90 {-.122008467928146216, -.211324865405187118 , .577350269189625764},
91 { .122008467928146216, -.211324865405187118 , .577350269189625764}};
92
93 static double wt4[] = {0.394337567297406440, 0.394337567297406440,
94 0.105662432702593559, 0.105662432702593559};
95
96 static double rst9[][3] = {
97 {mpt68, mpt88, mpt77},
98 { zero, mpt88, mpt77},
99 { pt68, mpt88, mpt77},
100 {mpt38, mpt5, zero},
101 { zero, mpt5, zero},
102 { pt38, mpt5, zero},
103 {mpt08, mpt11, pt77},
104 { zero, mpt11, pt77},
105 { pt08, mpt11, pt77}
106 };
107 static double wt9[] = {Qw311, Qw321, Qw331, Qw312, Qw322, Qw332, Qw313,
108 Qw323, Qw333};
109
110 /* static double rst7[][3] = {{A,A,A}, */
111 /* {O,P,P},{P,O,P},{P,P,O}, */
112 /* {Q,Q,R},{Q,R,Q},{R,Q,Q}}; */
113 /* static double wt7[] = {L,M,M,M,N,N,N}; */
114
115 /* static double rst12[][3] = {{V,W,1.0-V-W},{W,V,1.0-V-W},{W,W,1.0-W-W}, */
116 /* {X,Y,1.0-X-Y},{Y,X,1.0-X-Y},{Y,Y,1.0-Y-Y}, */
117 /* {Z,AA,1.0-Z-AA},{AA,Z,1.0-Z-AA},{Z,BB,1.0-Z-BB}, */
118 /* {AA,BB,1.0-AA-BB},{BB,AA,1.0-BB-AA},{BB,Z,1.0-BB-Z} */
119
120 /* }; */
121
122 /* static double wt12[] = {S,S,S,T,T,T,U,U,U,U,U,U}; */
123
triIntPntPyr(int nint,DARR3 ** bcord,double ** wt)124 int triIntPntPyr(int nint, DARR3 **bcord, double **wt)
125 {
126 int retval = 1 ;
127 if( nint == 1 ){*bcord = rst1 ; *wt = wt1; }
128 else if( nint == 4 ){*bcord = rst4 ; *wt = wt4; }
129 else if( nint == 9 ){*bcord = rst9 ; *wt = wt9; }
130 else
131 {
132 /* fprintf(stderr,"\n%d integration points unsupported in symtri.c; give {1,3,4,6,7,12}\n",nint); */
133 fprintf(stderr,"\n%d integration points unsupported in symtri.c; give {1,4}\n",nint);
134 retval = 0;
135 }
136 return retval ;
137 }
138
139
140
141