| 430 | } // constructPatchLevelSCLaplaceOp |
| 431 | |
| 432 | void |
| 433 | PETScMatUtilities::constructPatchLevelVCSCViscousOp( |
| 434 | Mat& mat, |
| 435 | const SAMRAI::solv::PoissonSpecifications& poisson_spec, |
| 436 | double alpha, |
| 437 | double beta, |
| 438 | const std::vector<SAMRAI::solv::RobinBcCoefStrategy<NDIM>*>& bc_coefs, |
| 439 | double data_time, |
| 440 | const std::vector<int>& num_dofs_per_proc, |
| 441 | int dof_index_idx, |
| 442 | SAMRAI::tbox::Pointer<SAMRAI::hier::PatchLevel<NDIM>> patch_level, |
| 443 | VCInterpType mu_interp_type) |
| 444 | { |
| 445 | #if !defined(NDEBUG) |
| 446 | TBOX_ASSERT(bc_coefs.size() == NDIM); |
| 447 | #endif |
| 448 | |
| 449 | int ierr; |
| 450 | if (mat) |
| 451 | { |
| 452 | ierr = MatDestroy(&mat); |
| 453 | IBTK_CHKERRQ(ierr); |
| 454 | } |
| 455 | |
| 456 | // Determine the index ranges. |
| 457 | const int mpi_rank = IBTK_MPI::getRank(); |
| 458 | const int n_local = num_dofs_per_proc[mpi_rank]; |
| 459 | const int proc_lower = std::accumulate(num_dofs_per_proc.begin(), num_dofs_per_proc.begin() + mpi_rank, 0); |
| 460 | const int proc_upper = proc_lower + n_local; |
| 461 | const int n_total = std::accumulate(num_dofs_per_proc.begin(), num_dofs_per_proc.end(), 0); |
| 462 | |
| 463 | // Determine the non-zero structure of the matrix. |
| 464 | std::vector<int> d_nnz(n_local, 0), o_nnz(n_local, 0); |
| 465 | for (PatchLevel<NDIM>::Iterator p(patch_level); p; p++) |
| 466 | { |
| 467 | Pointer<Patch<NDIM>> patch = patch_level->getPatch(p()); |
| 468 | const Box<NDIM>& patch_box = patch->getBox(); |
| 469 | Pointer<SideData<NDIM, int>> dof_index_data = patch->getPatchData(dof_index_idx); |
| 470 | #if !defined(NDEBUG) |
| 471 | TBOX_ASSERT(dof_index_data->getDepth() == 1); |
| 472 | #endif |
| 473 | for (unsigned int axis = 0; axis < NDIM; ++axis) |
| 474 | { |
| 475 | for (Box<NDIM>::Iterator b(SideGeometry<NDIM>::toSideBox(patch_box, axis)); b; b++) |
| 476 | { |
| 477 | const hier::Index<NDIM>& cc = b(); |
| 478 | const SideIndex<NDIM> i(cc, axis, SideIndex<NDIM>::Lower); |
| 479 | const int i_dof_index = (*dof_index_data)(i); |
| 480 | if (proc_lower <= i_dof_index && i_dof_index < proc_upper) |
| 481 | { |
| 482 | // Stencil for finite difference operator. |
| 483 | const int local_idx = i_dof_index - proc_lower; |
| 484 | d_nnz[local_idx] += 1; |
| 485 | |
| 486 | for (unsigned int d = 0; d < NDIM; ++d) |
| 487 | { |
| 488 | if (d == axis) |
| 489 | { |