AcouSTO  version 2.0

◆ intgba()

intgba ( int  kcrce,
int  kelem,
int  kode,
struct vector *  xcntr,
struct panel4 *  element,
DOUBLE *  source,
DOUBLE *  doublet,
DOUBLE *  ratelet 
)

Analytically evaluates quadrature on a quadrangular patch

Parameters
[in]kcrce
[in]kelem
[in]kode
[in]xcntrControl point
[in]elementPanel
[out]sourceSource coefficient
[out]doubletDoublet coefficient
[out]rateletRatelet coefficient.
48  {
49 
50  int i,l,m,n,ks;
51  int icsi,ieta;
52  int icornr,jcornr;
53  int iloop,nloop;
54  int korn[3][4]={{0,1,2,3},{0,1,1,3},{3,1,2,3}};
55  int kcsi[4]={0,0,1,1};
56  int keta[4]={1,0,0,1};
57  DOUBLE signc,signj[4]={-1.0,1.0,-1.0,1.0};
58  DOUBLE pi;
59  DOUBLE consta;
60  DOUBLE signe;
61  DOUBLE eps;
62  struct vector pij[4];
63  struct vector a1[2],a2[2];
64  struct vector a1p[2],a2p[2];
65  struct vector vtmp;
66  int natri;
67  DOUBLE a1a1 ,a2a2;
68  DOUBLE side1,side2;
69  struct vector unc,un[4];
70  DOUBLE unac;
71  DOUBLE bij,cij;
72  struct vector pz;
73  DOUBLE pzun,rc2,rc1;
74 
75  struct vector q[4];
76  struct vector qxa1;
77  struct vector qxa2;
78  DOUBLE qa1xa2[4];
79  DOUBLE sj;
80  COMPLEX cplt;
81  DOUBLE tanp,hnumer,qdotq,denom;
82  DOUBLE aa1,qdota1,aqxa1,alog1;
83  DOUBLE aa2,qdota2,aqxa2,alog2;
84  DOUBLE unqxa1,unqxa2;
85 
86  pi = 4*atan(1.0);
87  consta = .25/pi;
88  signe = (DOUBLE)kcrce;
89  eps = 0.0E0;
90 
91  /* Centroids and base vectors */
92  pij[0].x = .25*(element->x[0].x+ element->x[1].x+ element->x[2].x+ element->x[3].x);
93  pij[0].y = .25*(element->x[0].y+ element->x[1].y+ element->x[2].y+ element->x[3].y);
94  pij[0].z = .25*(element->x[0].z+ element->x[1].z+ element->x[2].z+ element->x[3].z);
95 
96  pij[1].x = .25*(element->x[1].x+ element->x[0].x- element->x[2].x- element->x[3].x);
97  pij[1].y = .25*(element->x[1].y+ element->x[0].y- element->x[2].y- element->x[3].y);
98  pij[1].z = .25*(element->x[1].z+ element->x[0].z- element->x[2].z- element->x[3].z);
99 
100  pij[2].x = .25*(element->x[1].x- element->x[0].x+ element->x[2].x- element->x[3].x);
101  pij[2].y = .25*(element->x[1].y- element->x[0].y+ element->x[2].y- element->x[3].y);
102  pij[2].z = .25*(element->x[1].z- element->x[0].z+ element->x[2].z- element->x[3].z);
103 
104  pij[3].x = .25*(element->x[1].x- element->x[0].x- element->x[2].x+ element->x[3].x);
105  pij[3].y = .25*(element->x[1].y- element->x[0].y- element->x[2].y+ element->x[3].y);
106  pij[3].z = .25*(element->x[1].z- element->x[0].z- element->x[2].z+ element->x[3].z);
107 
108 
109  vec_diff (&a1[0],element->x[1],element->x[2]); vec_scale(&a1[0],.5);
110  vec_diff (&a1[1],element->x[0],element->x[3]); vec_scale(&a1[1],.5);
111 
112  vec_diff (&a2[0],element->x[1],element->x[0]); vec_scale(&a2[0],.5);
113  vec_diff (&a2[1],element->x[2],element->x[3]); vec_scale(&a2[1],.5);
114 
115 
116  natri = 0;
117  side1 = 0.0;
118  side2 = 0.0;
119  /*
120  l 1 2 -> 0 1
121  m 2 1 -> 1 0
122  */
123  for(l=0;l<=1;l++){
124  m=1-l;
125  a1a1 = vec_mod_square(a1[l]);
126  side1 += sqrt(a1a1);
127  if( a1a1 > eps) goto l73;
128  natri +=1;
129  vec_copy(&a1[l],a1[m]);
130  l73:
131  a2a2 = vec_mod_square(a2[l]);
132  side2 += sqrt(a2a2);
133  if(a2a2 <= eps){
134  natri += 1;
135  vec_copy(&a2[l],a2[m]);
136  }
137  }
138 
139  /* unit normal */
140  unc.x = ((a1[0].y + a1[1].y)*(a2[0].z + a2[1].z) - (a1[0].z + a1[1].z)*(a2[0].y + a2[1].y))*.25;
141  unc.y = ((a1[0].z + a1[1].z)*(a2[0].x + a2[1].x) - (a1[0].x + a1[1].x)*(a2[0].z + a2[1].z))*.25;
142  unc.z = ((a1[0].x + a1[1].x)*(a2[0].y + a2[1].y) - (a1[0].y + a1[1].y)*(a2[0].x + a2[1].x))*.25;
143  unac = vec_mod(unc);
144  vec_normalize(&unc);
145 
146 
147  bij=0.0;
148  cij=0.0;
149  vec_diff(&pz, pij[0],*xcntr);
150  pzun = vec_dot(pz,unc);
151  rc2 = vec_mod_square(pz);
152  rc1 = sqrt(rc2);
153 
154  /* analytical calculation of source and doublet integrals */
155  for(i=0;i<4;i++){
156  vec_diff(&q[i],element->x[i],*xcntr);
157  }
158 
159  /*
160  iloop 1 2 3 0 1 2
161  m 0 1 2 -1 0 1
162  n 3 2 1 2 1 0
163  */
164  for(iloop=0;iloop<3;iloop++){
165  for(l=0;l<2;l++){
166  vec_copy(&a1p[l],a1[l]);
167  vec_copy(&a2p[l],a2[l]);
168  }
169  if(0 == iloop) goto l105;
170  m = iloop-1;
171  n = 2 -iloop;
172  vec_copy(&a1p[m],a1p[n]);
173  vec_zero(&vtmp);
174  vec_diff(&vtmp,element->x[1],element->x[3]);
175  vec_scale(&vtmp,.5);
176  vec_copy(&a2p[n], vtmp);
177 
178  l105:
179  ks=0;
180 
181  for(jcornr=0;jcornr<4;jcornr++){
182  icsi = kcsi[jcornr];
183  ieta = keta[jcornr];
184  icornr = korn[iloop][jcornr];
185  qa1xa2[jcornr] = 0.0;
186 
187  vec_cross(&un[jcornr],a1p[ieta],a2p[icsi]);
188  qa1xa2[jcornr] = vec_dot(q[icornr],un[jcornr]);
189 
190  if(qa1xa2[jcornr] < (-eps)){
191  ks= ks -1;
192  }
193  if(qa1xa2[jcornr] > (eps)){
194  ks+=1;
195  }
196  }
197 
198  nloop = 0;
199  if(0 != iloop ) goto l125;
200  nloop = 2;
201  if(1 == kode ) goto l125;
202  if(4 == abs(ks) || 0 != natri ) goto l125;
203  nloop = 0;
204  goto l155;
205 
206  l125:
207  sj = 0.0;
208  if(fabs(qa1xa2[0]) > eps) sj = SIGN(1,qa1xa2[1]);
209  for(jcornr=0;jcornr<4;jcornr++){
210  icsi = kcsi[jcornr];
211  ieta = keta[jcornr];
212  signc = signj[jcornr];
213  icornr = korn[iloop][jcornr];
214  vec_cross(&qxa1,q[icornr],a1p[ieta]);
215  vec_cross(&qxa2,q[icornr],a2p[icsi]);
216  tanp = 0.0;
217  if(4 != abs(ks) ) goto l140;
218  if(1 == kode ) goto l140;
219 
220  hnumer = vec_dot(qxa1,qxa2);
221  qdotq = vec_dot(q[icornr],q[icornr]);
222  denom = sqrt(qdotq) * fabs(qa1xa2[jcornr]);
223  cplt = hnumer + I*denom;
224  if(denom > eps) tanp = atan2l(hnumer,denom)*sj;
225  cij+=tanp * signc;
226  l140:
227 
228 
229  aa1 = vec_mod(a1p[ieta]);
230  qdota1 = vec_dot(q[icornr],a1p[ieta]);
231  aqxa1 = vec_mod(qxa1);
232  alog1 = 0.0;
233  if(aqxa1 > eps) alog1 = asinhl(qdota1/aqxa1)/aa1;
234 
235  aa2 = vec_mod(a2p[icsi]);
236  qdota2 = vec_dot(q[icornr],a2p[icsi]);
237  aqxa2 = vec_mod(qxa2);
238  alog2 = 0.0;
239  if(aqxa2 > eps) alog2 = asinhl(qdota2/aqxa2)/aa2;
240 
241  unqxa1 = vec_dot(unc,qxa1);
242  unqxa2 = vec_dot(unc,qxa2);
243  bij = bij - signc*(-unqxa1*alog1+unqxa2*alog2+pzun*tanp);
244 
245  }
246 
247  l155:
248  if(2 == nloop) break;
249  }
250 
251  *source = bij*consta;
252  *doublet = cij*consta*signe;
253  *ratelet = cij*consta*signe*rc1;
254 
255 
256 }
vec_normalize
void vec_normalize(struct vector *v)
Definition: math.c:130
vec_scale
void vec_scale(struct vector *v1, DOUBLE d)
Definition: math.c:170
vec_copy
void vec_copy(struct vector *vdest, const struct vector vsrc)
Definition: math.c:158
SIGN
#define SIGN(x, y)
Definition: functions.h:27
vector
vector struct to hold triplets.
Definition: structs.h:29
vec_mod
DOUBLE vec_mod(struct vector v1)
Definition: math.c:121
vec_zero
void vec_zero(struct vector *v)
Definition: math.c:146
COMPLEX
#define COMPLEX
Definition: types.h:48
vec_mod_square
DOUBLE vec_mod_square(struct vector v1)
Definition: math.c:108
vector::z
DOUBLE z
Definition: structs.h:35
vec_diff
void vec_diff(struct vector *vdest, const struct vector v1, const struct vector v2)
Definition: math.c:52
vec_dot
DOUBLE vec_dot(struct vector v1, struct vector v2)
Definition: math.c:34
DOUBLE
double DOUBLE
Definition: types.h:44
vec_cross
void vec_cross(struct vector *vdest, const struct vector v1, const struct vector v2)
Definition: math.c:92
vector::x
DOUBLE x
Definition: structs.h:31
panel4::x
struct vector x[4]
Definition: structs.h:46
vector::y
DOUBLE y
Definition: structs.h:33