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 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 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