219 for (
int m = 0; m < 9; ++m) {
223 auto& mat_tmp =
reinterpret_cast<Array2D<Real,-1,1,-1,1
>&>(sten);
227 if (flag(i,j,k).isRegular())
229 mat_tmp(0,0) = fac[0]*(
b[0](i,j,k)+
b[0](i+1,j,k))
230 + fac[1]*(
b[1](i,j,k)+
b[1](i,j+1,k));
231 mat_tmp(-1, 0) = -fac[0]*
b[0](i,j,k);
232 mat_tmp( 1, 0) = -fac[0]*
b[0](i+1,j,k);
233 mat_tmp( 0,-1) = -fac[1]*
b[1](i,j,k);
234 mat_tmp( 0, 1) = -fac[1]*
b[1](i,j+1,k);
236 if (cell_id(i-1,j,k) < 0) {
240 mat_tmp(0,0) += bf1*
b[0](i,j,k);
241 mat_tmp(-1,0) = Real(0.0);
242 mat_tmp(1,0) += bf2*
b[0](i,j,k);
245 if (cell_id(i+1,j,k) < 0) {
249 mat_tmp(0,0) += bf1*
b[0](i+1,j,k);
250 mat_tmp(1,0) = Real(0.0);
251 mat_tmp(-1,0) += bf2*
b[0](i+1,j,k);
254 if (cell_id(i,j-1,k) < 0) {
258 mat_tmp(0,0) += bf1*
b[1](i,j,k);
259 mat_tmp(0,-1) = Real(0.0);
260 mat_tmp(0,1) += bf2*
b[1](i,j,k);
263 if (cell_id(i,j+1,k) < 0) {
267 mat_tmp(0,0) += bf1*
b[1](i,j+1,k);
268 mat_tmp(0,1) = Real(0.0);
269 mat_tmp(0,-1) += bf2*
b[1](i,j+1,k);
272 if (sa != Real(0.0)) {
273 mat_tmp(0,0) += sa*a(i,j,k);
276 else if (flag(i,j,k).isSingleValued())
279 Real area =
apx(i,j,k);
281 if (area > Real(0.0))
287 if (area != Real(1.0)) {
288 joff =
static_cast<int>(std::copysign(Real(1.0),
fcx(i,j,k)));
290 if (cell_id(i-1,jj,k) < 0 && cell_id(i,jj,k) < 0) {
293 fracy = std::abs(
fcx(i,j,k));
304 Real tmp = (Real(1.0)-fracy)*area*
b[0](i,j,k);
305 if (cell_id(i-1,j,k) >= 0) {
306 mat_tmp( 0,0) += tmp*f;
307 mat_tmp(-1,0) -= tmp*f;
308 }
else if (cell_id(i+1,j,k) < 0 ||
apx(i+1,j,k) == Real(0.0)) {
309 mat_tmp(0,0) += tmp*(f+bflo);
311 mat_tmp(0,0) += tmp*(f+bf1);
312 mat_tmp(1,0) += tmp* bf2;
315 if (fracy > Real(0.0)) {
316 if (cell_id(i-1,jj,k) >= 0 && cell_id(i,jj,k) >= 0) {
317 mat_tmp(-1,joff) -= fracy*area*
b[0](i,jj,k)*f;
318 mat_tmp( 0,joff) += fracy*area*
b[0](i,jj,k)*f;
319 }
else if (cell_id(i+1,jj,k) < 0 ||
apx(i+1,jj,k) == Real(0.0)) {
320 mat_tmp( 0,joff) += tmp*(f+bflo);
322 mat_tmp( 0,joff) += tmp*(f+bf1);
323 mat_tmp( 1,joff) += tmp* bf2;
330 if (area > Real(0.0))
336 if (area != Real(1.0)) {
337 joff =
static_cast<int>(std::copysign(Real(1.0),
fcx(i+1,j,k)));
339 if (cell_id(i,jj,k) < 0 && cell_id(i+1,jj,k) < 0) {
342 fracy = std::abs(
fcx(i+1,j,k));
353 Real tmp = (Real(1.0)-fracy)*area*
b[0](i+1,j,k);
354 if (cell_id(i+1,j,k) >= 0) {
355 mat_tmp(0,0) += tmp*f;
356 mat_tmp(1,0) -= tmp*f;
357 }
else if (cell_id(i-1,j,k) < 0 ||
apx(i,j,k) == Real(0.0)) {
358 mat_tmp(0,0) += tmp*(f+bflo);
360 mat_tmp( 0,0) += tmp*(f+bf1);
361 mat_tmp(-1,0) += tmp* bf2;
364 if (fracy > Real(0.0)) {
365 if (cell_id(i,jj,k) >= 0 && cell_id(i+1,jj,k) >= 0) {
366 mat_tmp(0,joff) += fracy*area*
b[0](i+1,jj,k)*f;
367 mat_tmp(1,joff) -= fracy*area*
b[0](i+1,jj,k)*f;
368 }
else if (cell_id(i-1,jj,k) < 0 ||
apx(i,jj,k) == Real(0.0)) {
369 mat_tmp(0,joff) += tmp*(f+bflo);
371 mat_tmp( 0,joff) += tmp*(f+bf1);
372 mat_tmp(-1,joff) += tmp* bf2;
379 if (area > Real(0.0))
385 if (area != Real(1.0)) {
386 ioff =
static_cast<int>(std::copysign(Real(1.0),
fcy(i,j,k)));
388 if (cell_id(ii,j-1,k) < 0 && cell_id(ii,j,k) < 0) {
391 fracx = std::abs(
fcy(i,j,k));
402 Real tmp = (Real(1.0)-fracx)*area*
b[1](i,j,k);
403 if (cell_id(i,j-1,k) >= 0) {
404 mat_tmp(0, 0) += tmp*f;
405 mat_tmp(0,-1) -= tmp*f;
406 }
else if (cell_id(i,j+1,k) < 0 ||
apy(i,j+1,k) == Real(0.0)) {
407 mat_tmp(0,0) += tmp*(f+bflo);
409 mat_tmp(0,0) += tmp*(f+bf1);
410 mat_tmp(0,1) += tmp* bf2;
413 if (fracx > Real(0.0)) {
414 if (cell_id(ii,j-1,k) >= 0 && cell_id(ii,j,k) >= 0) {
415 mat_tmp(ioff,-1) -= fracx*area*
b[1](ii,j,k)*f;
416 mat_tmp(ioff, 0) += fracx*area*
b[1](ii,j,k)*f;
417 }
else if (cell_id(ii,j+1,k) < 0 ||
apy(ii,j+1,k) == Real(0.0)) {
418 mat_tmp(ioff,0) += tmp*(f+bflo);
420 mat_tmp(ioff,0) += tmp*(f+bf1);
421 mat_tmp(ioff,1) += tmp* bf2;
428 if (area > Real(0.0))
434 if (area != Real(1.0)) {
435 ioff =
static_cast<int>(std::copysign(Real(1.0),
fcy(i,j+1,k)));
437 if (cell_id(ii,j,k) < 0 && cell_id(ii,j+1,k) < 0) {
440 fracx = std::abs(
fcy(i,j+1,k));
451 Real tmp = (Real(1.0)-fracx)*area*
b[1](i,j+1,k);
452 if (cell_id(i,j+1,k) >= 0) {
453 mat_tmp(0,0) += tmp*f;
454 mat_tmp(0,1) -= tmp*f;
455 }
else if (cell_id(i,j-1,k) < 0 ||
apy(i,j,k) == Real(0.0)) {
456 mat_tmp(0,0) += tmp*(f+bflo);
458 mat_tmp(0, 0) += tmp*(f+bf1);
459 mat_tmp(0,-1) += tmp* bf2;
462 if (fracx > Real(0.0)) {
463 if (cell_id(ii,j,k) >= 0 && cell_id(ii,j+1,k) >= 0) {
464 mat_tmp(ioff,0) += fracx*area*
b[1](ii,j+1,k)*f;
465 mat_tmp(ioff,1) -= fracx*area*
b[1](ii,j+1,k)*f;
466 }
else if (cell_id(ii,j-1,k) < 0 ||
apy(ii,j,k) == Real(0.0)) {
467 mat_tmp(ioff,0) += tmp*(f+bflo);
469 mat_tmp(ioff, 0) += tmp*(f+bf1);
470 mat_tmp(ioff,-1) += tmp* bf2;
476 Real anorm = std::sqrt((
apx(i,j,k) -
apx(i+1,j,k))*(
apx(i,j,k) -
apx(i+1,j,k))
477 + (
apy(i,j,k) -
apy(i,j+1,k))*(
apy(i,j,k) -
apy(i,j+1,k)));
479 Real anorminv = Real(1.0)/anorm;
480 Real anrmx = (
apx(i,j,k) -
apx(i+1,j,k))*anorminv;
481 Real anrmy = (
apy(i,j,k) -
apy(i,j+1,k))*anorminv;
482 Real sx = std::copysign(Real(1.0),anrmx);
483 Real sy = std::copysign(Real(1.0),anrmy);
484 Real bctx = bcen(i,j,k,0);
485 Real bcty = bcen(i,j,k,1);
487 Real gx = bctx - dg*anrmx;
488 Real gy = bcty - dg*anrmy;
489 int ioff = -
static_cast<int>(sx);
490 int joff = -
static_cast<int>(sy);
492 - gx*sx - gx*gy*sx*sy,
493 - gy*sy - gx*gy*sx*sy,
497 for (
int ii=0; ii<4; ii++){
498 feb(ii) = -phig1(ii) * (ba(i,j,k) * beb(i,j,k) / dg);
500 mat_tmp(0 , 0 ) -= feb(0)*fac[0];
501 mat_tmp(ioff, 0 ) -= feb(1)*fac[0];
502 mat_tmp(0 ,joff) -= feb(2)*fac[0];
503 mat_tmp(ioff,joff) -= feb(3)*fac[0];
506 for (
int jj=-1; jj<=1; jj++) {
507 for (
int ii=-1; ii<=1; ii++) {
508 mat_tmp(ii,jj) *= Real(1.0)/vfrc(i,j,k);
511 if (sa != Real(0.0)) {
512 mat_tmp(0,0) += sa*a(i,j,k);
517 for (
int jj=-1; jj<=1; jj++) {
518 for (
int ii=-1; ii<=1; ii++) {
519 mat_tmp(ii,jj) = Real(0.0);
521 mat_tmp(0,0) = Real(1.0);
524 diaginv(i,j,k) = Real(1.0) / mat_tmp(0,0);
525 for (
int jj=-1; jj<=1; jj++) {
526 for (
int ii=-1; ii<=1; ii++) {
527 mat_tmp(ii,jj) *= diaginv(i,j,k);
529 mat_tmp(0,0) = Real(1.0);
532 for (
int m = 0; m < 9; ++m) {
533 ncols(i,j,k) += (sten[m] != Real(0.0));