| 126 |
|
return(sgn2 == sgn3); |
| 127 |
|
} |
| 128 |
|
|
| 129 |
< |
/* Test and set for edge */ |
| 129 |
> |
/* Test (and set) bitmap for edge */ |
| 130 |
|
static int |
| 131 |
|
check_edge(unsigned char *emap, int nedges, const MIGRATION *mig, int mark) |
| 132 |
|
{ |
| 202 |
|
input_orient*invec[2]) { |
| 203 |
|
for (miga[0] = rbf->ejl; miga[0] != NULL; |
| 204 |
|
miga[0] = nextedge(rbf,miga[0])) |
| 205 |
< |
if (opp_rbf(rbf,miga[0]) == rbf->next) |
| 205 |
> |
if (opp_rbf(rbf,miga[0]) == rbf->next) { |
| 206 |
> |
double nf = 1.-rbf->invec[2]*rbf->invec[2]; |
| 207 |
> |
if (nf > FTINY) { |
| 208 |
> |
nf = sqrt((1.-invec[2]*invec[2])/nf); |
| 209 |
> |
invec[0] = nf*rbf->invec[0]; |
| 210 |
> |
invec[1] = nf*rbf->invec[1]; |
| 211 |
> |
} |
| 212 |
|
return(0); |
| 213 |
+ |
} |
| 214 |
|
break; |
| 215 |
|
} |
| 216 |
|
} |
| 254 |
|
double t, full_dist; |
| 255 |
|
/* get relative position */ |
| 256 |
|
t = acos(DOT(invec, mig->rbfv[0]->invec)); |
| 257 |
< |
if (t < M_PI/GRIDRES) { /* near first DSF */ |
| 257 |
> |
if (t < M_PI/grid_res) { /* near first DSF */ |
| 258 |
|
n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1); |
| 259 |
|
rbf = (RBFNODE *)malloc(n); |
| 260 |
|
if (rbf == NULL) |
| 263 |
|
return(rbf); |
| 264 |
|
} |
| 265 |
|
full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec)); |
| 266 |
< |
if (t > full_dist-M_PI/GRIDRES) { /* near second DSF */ |
| 266 |
> |
if (t > full_dist-M_PI/grid_res) { /* near second DSF */ |
| 267 |
|
n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1); |
| 268 |
|
rbf = (RBFNODE *)malloc(n); |
| 269 |
|
if (rbf == NULL) |
| 305 |
|
rbf->rbfa[n].crad = ANG2R(sqrt(rad0*rad0*(1.-t) + |
| 306 |
|
rad1*rad1*t)); |
| 307 |
|
ovec_from_pos(v, rbf1j->gx, rbf1j->gy); |
| 308 |
< |
geodesic(v, v0, v, t, GEOD_REL); |
| 308 |
> |
geodesic(v, v0, v, t*full_dist, GEOD_RAD); |
| 309 |
|
pos_from_vec(pos, v); |
| 310 |
|
rbf->rbfa[n].gx = pos[0]; |
| 311 |
|
rbf->rbfa[n].gy = pos[1]; |
| 331 |
|
float mbfact, mcfact; |
| 332 |
|
int n, i, j, k; |
| 333 |
|
FVECT v0, v1, v2; |
| 334 |
< |
double s, t; |
| 334 |
> |
double s, t, s_full, t_full; |
| 335 |
|
|
| 336 |
|
VCOPY(sivec, invec); /* find triangle/edge */ |
| 337 |
|
sym = get_interp(miga, sivec); |
| 339 |
|
return(NULL); |
| 340 |
|
if (miga[1] == NULL) { /* advect along edge? */ |
| 341 |
|
rbf = e_advect_rbf(miga[0], sivec); |
| 342 |
< |
rev_rbf_symmetry(rbf, sym); |
| 342 |
> |
if (single_plane_incident) |
| 343 |
> |
rotate_rbf(rbf, invec); |
| 344 |
> |
else |
| 345 |
> |
rev_rbf_symmetry(rbf, sym); |
| 346 |
|
return(rbf); |
| 347 |
|
} |
| 348 |
|
#ifdef DEBUG |
| 360 |
|
normalize(v2); |
| 361 |
|
fcross(v1, sivec, miga[1]->rbfv[1]->invec); |
| 362 |
|
normalize(v1); |
| 363 |
< |
s = acos(DOT(v0,v1)) / acos(DOT(v0,v2)); |
| 363 |
> |
s = acos(DOT(v0,v1)); |
| 364 |
|
geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec, |
| 365 |
< |
s, GEOD_REL); |
| 366 |
< |
t = acos(DOT(v1,sivec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec)); |
| 365 |
> |
s, GEOD_RAD); |
| 366 |
> |
s /= s_full = acos(DOT(v0,v2)); |
| 367 |
> |
t = acos(DOT(v1,sivec)) / |
| 368 |
> |
(t_full = acos(DOT(v1,miga[1]->rbfv[1]->invec))); |
| 369 |
|
n = 0; /* count migrating particles */ |
| 370 |
|
for (i = 0; i < mtx_nrows(miga[0]); i++) |
| 371 |
|
for (j = 0; j < mtx_ncols(miga[0]); j++) |
| 406 |
|
rad1j = R2ANG(rbf1j->crad); |
| 407 |
|
srad2 = (1.-s)*(1.-t)*rad0i*rad0i + s*(1.-t)*rad1j*rad1j; |
| 408 |
|
ovec_from_pos(v1, rbf1j->gx, rbf1j->gy); |
| 409 |
< |
geodesic(v1, v0, v1, s, GEOD_REL); |
| 409 |
> |
geodesic(v1, v0, v1, s*s_full, GEOD_RAD); |
| 410 |
|
for (k = 0; k < mtx_ncols(miga[2]); k++) { |
| 411 |
|
float mb = mtx_coef(miga[1],j,k); |
| 412 |
|
float mc = mtx_coef(miga[2],i,k); |
| 421 |
|
rad2k = R2ANG(rbf2k->crad); |
| 422 |
|
rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k)); |
| 423 |
|
ovec_from_pos(v2, rbf2k->gx, rbf2k->gy); |
| 424 |
< |
geodesic(vout, v1, v2, t, GEOD_REL); |
| 424 |
> |
geodesic(vout, v1, v2, t*t_full, GEOD_RAD); |
| 425 |
|
pos_from_vec(pos, vout); |
| 426 |
|
rbf->rbfa[n].gx = pos[0]; |
| 427 |
|
rbf->rbfa[n].gy = pos[1]; |