自学教程：C++ IBTK_CHKERRQ函数代码示例

voidPenaltyIBMethod::preprocessIntegrateData(double current_time, double new_time, int num_cycles){    IBMethod::preprocessIntegrateData(current_time, new_time, num_cycles);    const int coarsest_ln = 0;    const int finest_ln = d_hierarchy->getFinestLevelNumber();    // Look-up or allocate Lagangian data.    d_K_data.resize(finest_ln + 1);    d_M_data.resize(finest_ln + 1);    d_Y_current_data.resize(finest_ln + 1);    d_Y_new_data.resize(finest_ln + 1);    d_V_current_data.resize(finest_ln + 1);    d_V_new_data.resize(finest_ln + 1);    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)    {        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;        d_K_data[ln] = d_l_data_manager->getLData("K", ln);        d_M_data[ln] = d_l_data_manager->getLData("M", ln);        d_Y_current_data[ln] = d_l_data_manager->getLData("Y", ln);        d_Y_new_data[ln] = d_l_data_manager->createLData("Y_new", ln, NDIM);        d_V_current_data[ln] = d_l_data_manager->getLData("V", ln);        d_V_new_data[ln] = d_l_data_manager->createLData("V_new", ln, NDIM);        // Initialize Y^{n+1} and V^{n+1} to equal Y^{n} and V^{n}.        int ierr;        ierr = VecCopy(d_Y_current_data[ln]->getVec(), d_Y_new_data[ln]->getVec());        IBTK_CHKERRQ(ierr);        ierr = VecCopy(d_V_current_data[ln]->getVec(), d_V_new_data[ln]->getVec());        IBTK_CHKERRQ(ierr);    }    return;} // preprocessIntegrateData

voidIBImplicitModHelmholtzPETScLevelSolver::deallocateSolverState(){    if (!d_is_initialized) return;    IBAMR_TIMER_START(t_deallocate_solver_state);    // Deallocate PETSc objects.    int ierr;    ierr = KSPDestroy(d_petsc_ksp); IBTK_CHKERRQ(ierr);    ierr = MatDestroy(d_petsc_mat); IBTK_CHKERRQ(ierr);    ierr = VecDestroy(d_petsc_x); IBTK_CHKERRQ(ierr);    ierr = VecDestroy(d_petsc_b); IBTK_CHKERRQ(ierr);    d_dof_index_fill.setNull();    d_petsc_ksp = PETSC_NULL;    d_petsc_mat = PETSC_NULL;    d_petsc_x = PETSC_NULL;    d_petsc_b = PETSC_NULL;    // Deallocate DOF index data.    Pointer<PatchLevel<NDIM> > level = d_hierarchy->getPatchLevel(d_level_num);    if (level->checkAllocated(d_dof_index_idx)) level->deallocatePatchData(d_dof_index_idx);    // Indicate that the solver is NOT initialized.    d_is_initialized = false;    IBAMR_TIMER_STOP(t_deallocate_solver_state);    return;}// deallocateSolverState

voidPenaltyIBMethod::postprocessIntegrateData(double current_time, double new_time, int num_cycles){    IBMethod::postprocessIntegrateData(current_time, new_time, num_cycles);    const int coarsest_ln = 0;    const int finest_ln = d_hierarchy->getFinestLevelNumber();    // Reset time-dependent Lagrangian data.    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)    {        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;        int ierr;        ierr = VecSwap(d_Y_current_data[ln]->getVec(), d_Y_new_data[ln]->getVec());        IBTK_CHKERRQ(ierr);        ierr = VecSwap(d_V_current_data[ln]->getVec(), d_V_new_data[ln]->getVec());        IBTK_CHKERRQ(ierr);    }    // Deallocate Lagrangian scratch data.    d_K_data.clear();    d_M_data.clear();    d_Y_current_data.clear();    d_Y_new_data.clear();    d_V_current_data.clear();    d_V_new_data.clear();    return;} // postprocessIntegrateData

void CCPoissonPETScLevelSolver::initializeSolverStateSpecialized(const SAMRAIVectorReal<NDIM, double>& x,                                                                 const SAMRAIVectorReal<NDIM, double>& /*b*/){    // Allocate DOF index data.    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();    const int x_idx = x.getComponentDescriptorIndex(0);    Pointer<CellDataFactory<NDIM, double> > x_fac = var_db->getPatchDescriptor()->getPatchDataFactory(x_idx);    const int depth = x_fac->getDefaultDepth();    Pointer<CellDataFactory<NDIM, int> > dof_index_fac =        var_db->getPatchDescriptor()->getPatchDataFactory(d_dof_index_idx);    dof_index_fac->setDefaultDepth(depth);    Pointer<PatchLevel<NDIM> > level = d_hierarchy->getPatchLevel(d_level_num);    if (!level->checkAllocated(d_dof_index_idx)) level->allocatePatchData(d_dof_index_idx);    // Setup PETSc objects.    int ierr;    PETScVecUtilities::constructPatchLevelDOFIndices(d_num_dofs_per_proc, d_dof_index_idx, level);    const int mpi_rank = SAMRAI_MPI::getRank();    ierr = VecCreateMPI(PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_x);    IBTK_CHKERRQ(ierr);    ierr = VecCreateMPI(PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_b);    IBTK_CHKERRQ(ierr);    PETScMatUtilities::constructPatchLevelCCLaplaceOp(        d_petsc_mat, d_poisson_spec, d_bc_coefs, d_solution_time, d_num_dofs_per_proc, d_dof_index_idx, level);    d_petsc_pc = d_petsc_mat;    d_petsc_ksp_ops_flag = SAME_PRECONDITIONER;    d_data_synch_sched = PETScVecUtilities::constructDataSynchSchedule(x_idx, level);    d_ghost_fill_sched = PETScVecUtilities::constructGhostFillSchedule(x_idx, level);    return;} // initializeSolverStateSpecialized

PetscErrorCode IBImplicitStaggeredHierarchyIntegrator::compositeIBJacobianApply(Vec x, Vec f){    PetscErrorCode ierr;    const double half_time = d_integrator_time + 0.5 * d_current_dt;    Vec* component_sol_vecs;    Vec* component_rhs_vecs;    ierr = VecMultiVecGetSubVecs(x, &component_sol_vecs);    IBTK_CHKERRQ(ierr);    ierr = VecMultiVecGetSubVecs(f, &component_rhs_vecs);    IBTK_CHKERRQ(ierr);    Pointer<SAMRAIVectorReal<NDIM, double> > u =        PETScSAMRAIVectorReal::getSAMRAIVector(component_sol_vecs[0]);    Pointer<SAMRAIVectorReal<NDIM, double> > f_u =        PETScSAMRAIVectorReal::getSAMRAIVector(component_rhs_vecs[0]);    Pointer<Variable<NDIM> > u_var = d_ins_hier_integrator->getVelocityVariable();    const int u_idx = u->getComponentDescriptorIndex(0);    const int f_u_idx = f_u->getComponentDescriptorIndex(0);    Vec X = component_sol_vecs[1];    Vec R = component_rhs_vecs[1];    // Evaluate the Eulerian terms.    d_stokes_op->setHomogeneousBc(true);    d_stokes_op->apply(*u, *f_u);    d_ib_implicit_ops->computeLinearizedLagrangianForce(X, half_time);    if (d_enable_logging)        plog << d_object_name             << "::integrateHierarchy(): spreading Lagrangian force to the Eulerian grid/n";    d_hier_velocity_data_ops->setToScalar(d_f_idx, 0.0);    d_u_phys_bdry_op->setPatchDataIndex(d_f_idx);    d_ib_implicit_ops->spreadLinearizedForce(d_f_idx,                                             d_u_phys_bdry_op,                                             getProlongRefineSchedules(d_object_name + "::f"),                                             half_time);    d_hier_velocity_data_ops->subtract(f_u_idx, f_u_idx, d_f_idx);    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(component_rhs_vecs[0]));    IBTK_CHKERRQ(ierr);    // Evaluate the Lagrangian terms.    d_hier_velocity_data_ops->scale(d_u_idx, 0.5, u_idx);    d_u_phys_bdry_op->setPatchDataIndex(d_u_idx);    d_ib_implicit_ops->interpolateLinearizedVelocity(        d_u_idx,        getCoarsenSchedules(d_object_name + "::u::CONSERVATIVE_COARSEN"),        getGhostfillRefineSchedules(d_object_name + "::u"),        half_time);    d_ib_implicit_ops->computeLinearizedResidual(X, R);    // Ensure that PETSc sees that the state of the RHS vector has changed.    // This is a nasty hack.    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(f));    IBTK_CHKERRQ(ierr);    return ierr;} // compositeIBJacobianApply

boolPETScKrylovLinearSolver::solveSystem(SAMRAIVectorReal<NDIM, double>& x, SAMRAIVectorReal<NDIM, double>& b){    IBTK_TIMER_START(t_solve_system);#if !defined(NDEBUG)    TBOX_ASSERT(d_A);#endif    int ierr;    // Initialize the solver, when necessary.    const bool deallocate_after_solve = !d_is_initialized;    if (deallocate_after_solve) initializeSolverState(x, b);#if !defined(NDEBUG)    TBOX_ASSERT(d_petsc_ksp);#endif    resetKSPOptions();    // Allocate scratch data.    d_b->allocateVectorData();    // Solve the system using a PETSc KSP object.    d_b->copyVector(Pointer<SAMRAIVectorReal<NDIM, double> >(&b, false));    d_A->setHomogeneousBc(d_homogeneous_bc);    d_A->modifyRhsForBcs(*d_b);    d_A->setHomogeneousBc(true);    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_x, Pointer<SAMRAIVectorReal<NDIM, double> >(&x, false));    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_b, d_b);    ierr = KSPSolve(d_petsc_ksp, d_petsc_b, d_petsc_x);    IBTK_CHKERRQ(ierr);    d_A->setHomogeneousBc(d_homogeneous_bc);    d_A->imposeSolBcs(x);    // Get iterations count and residual norm.    ierr = KSPGetIterationNumber(d_petsc_ksp, &d_current_iterations);    IBTK_CHKERRQ(ierr);    ierr = KSPGetResidualNorm(d_petsc_ksp, &d_current_residual_norm);    IBTK_CHKERRQ(ierr);    d_A->setHomogeneousBc(d_homogeneous_bc);    // Determine the convergence reason.    KSPConvergedReason reason;    ierr = KSPGetConvergedReason(d_petsc_ksp, &reason);    IBTK_CHKERRQ(ierr);    const bool converged = (static_cast<int>(reason) > 0);    if (d_enable_logging) reportKSPConvergedReason(reason, plog);    // Dealocate scratch data.    d_b->deallocateVectorData();    // Deallocate the solver, when necessary.    if (deallocate_after_solve) deallocateSolverState();    IBTK_TIMER_STOP(t_solve_system);    return converged;} // solveSystem

PetscErrorCodeIBImplicitStaggeredHierarchyIntegrator::compositeIBPCApply_SAMRAI(PC pc, Vec x, Vec y){    PetscErrorCode ierr;    void* ctx;    ierr = PCShellGetContext(pc, &ctx);    IBTK_CHKERRQ(ierr);    IBImplicitStaggeredHierarchyIntegrator* ib_integrator =        static_cast<IBImplicitStaggeredHierarchyIntegrator*>(ctx);    ierr = ib_integrator->compositeIBPCApply(x, y);    IBTK_CHKERRQ(ierr);    return ierr;} // compositeIBPCApply_SAMRAI

void GeneralizedIBMethod::preprocessIntegrateData(double current_time,                                                  double new_time,                                                  int num_cycles){    d_ib_force_and_torque_fcn_needs_init =        d_ib_force_fcn_needs_init || d_ib_force_and_torque_fcn_needs_init;    IBMethod::preprocessIntegrateData(current_time, new_time, num_cycles);    const int coarsest_ln = 0;    const int finest_ln = d_hierarchy->getFinestLevelNumber();    const double start_time = d_ib_solver->getStartTime();    if (d_ib_force_and_torque_fcn)    {        if (d_ib_force_and_torque_fcn_needs_init)        {            const bool initial_time =                MathUtilities<double>::equalEps(current_time, start_time);            resetLagrangianForceAndTorqueFunction(current_time, initial_time);            d_ib_force_and_torque_fcn_needs_init = false;        }    }    // Look-up or allocate Lagangian data.    d_D_current_data.resize(finest_ln + 1);    d_D_new_data.resize(finest_ln + 1);    d_N_current_data.resize(finest_ln + 1);    d_N_new_data.resize(finest_ln + 1);    d_W_current_data.resize(finest_ln + 1);    d_W_new_data.resize(finest_ln + 1);    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)    {        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;        d_D_current_data[ln] = d_l_data_manager->getLData("D", ln);        d_D_new_data[ln] = d_l_data_manager->createLData("D_new", ln, NDIM * NDIM);        d_N_current_data[ln] = d_l_data_manager->createLData("N", ln, NDIM);        d_N_new_data[ln] = d_l_data_manager->createLData("N_new", ln, NDIM);        d_W_current_data[ln] = d_l_data_manager->getLData("W", ln);        d_W_new_data[ln] = d_l_data_manager->createLData("W_new", ln, NDIM);        // Initialize D^{n+1} to equal D^{n}, and initialize W^{n+1} to equal        // W^{n}.        int ierr;        ierr = VecCopy(d_D_current_data[ln]->getVec(), d_D_new_data[ln]->getVec());        IBTK_CHKERRQ(ierr);        ierr = VecCopy(d_W_current_data[ln]->getVec(), d_W_new_data[ln]->getVec());        IBTK_CHKERRQ(ierr);    }    return;} // preprocessIntegrateData

void StaggeredStokesPETScLevelSolver::initializeSolverStateSpecialized(    const SAMRAIVectorReal<NDIM, double>& x,    const SAMRAIVectorReal<NDIM, double>& /*b*/){    // Allocate DOF index data.    Pointer<PatchLevel<NDIM> > level = d_hierarchy->getPatchLevel(d_level_num);    if (!level->checkAllocated(d_u_dof_index_idx)) level->allocatePatchData(d_u_dof_index_idx);    if (!level->checkAllocated(d_p_dof_index_idx)) level->allocatePatchData(d_p_dof_index_idx);    // Setup PETSc objects.    int ierr;    StaggeredStokesPETScVecUtilities::constructPatchLevelDOFIndices(        d_num_dofs_per_proc, d_u_dof_index_idx, d_p_dof_index_idx, level);    const int mpi_rank = SAMRAI_MPI::getRank();    ierr = VecCreateMPI(               PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_x);    IBTK_CHKERRQ(ierr);    ierr = VecCreateMPI(               PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_b);    IBTK_CHKERRQ(ierr);    StaggeredStokesPETScMatUtilities::constructPatchLevelMACStokesOp(d_petsc_mat,            d_U_problem_coefs,            d_U_bc_coefs,            d_new_time,            d_num_dofs_per_proc,            d_u_dof_index_idx,            d_p_dof_index_idx,            level);    ierr = MatDuplicate(d_petsc_mat, MAT_COPY_VALUES, &d_petsc_pc);    IBTK_CHKERRQ(ierr);    HierarchyDataOpsManager<NDIM>* hier_ops_manager =        HierarchyDataOpsManager<NDIM>::getManager();    Pointer<HierarchyDataOpsInteger<NDIM> > hier_p_dof_index_ops =        hier_ops_manager->getOperationsInteger(d_p_dof_index_var, d_hierarchy, true);    hier_p_dof_index_ops->resetLevels(d_level_num, d_level_num);    const int min_p_idx = hier_p_dof_index_ops->min(                              d_p_dof_index_idx); // NOTE: HierarchyDataOpsInteger::max() is broken    ierr = MatZeroRowsColumns(d_petsc_pc, 1, &min_p_idx, 1.0, NULL, NULL);    IBTK_CHKERRQ(ierr);    d_petsc_ksp_ops_flag = SAME_PRECONDITIONER;    const int u_idx = x.getComponentDescriptorIndex(0);    const int p_idx = x.getComponentDescriptorIndex(1);    d_data_synch_sched =        StaggeredStokesPETScVecUtilities::constructDataSynchSchedule(u_idx, p_idx, level);    d_ghost_fill_sched =        StaggeredStokesPETScVecUtilities::constructGhostFillSchedule(u_idx, p_idx, level);    return;} // initializeSolverStateSpecialized

PetscErrorCode VecAYPX_SAMRAI(Vec y, const PetscScalar alpha, Vec x){    IBTK_TIMER_START(t_vec_aypx);#if !defined(NDEBUG)    TBOX_ASSERT(x);    TBOX_ASSERT(y);#endif    static const bool interior_only = false;    if (MathUtilities<double>::equalEps(alpha, 1.0))    {        PSVR_CAST2(y)->add(PSVR_CAST2(x), PSVR_CAST2(y), interior_only);    }    else if (MathUtilities<double>::equalEps(alpha, -1.0))    {        PSVR_CAST2(y)->subtract(PSVR_CAST2(x), PSVR_CAST2(y), interior_only);    }    else    {        PSVR_CAST2(y)->axpy(alpha, PSVR_CAST2(y), PSVR_CAST2(x), interior_only);    }    int ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(y));    IBTK_CHKERRQ(ierr);    IBTK_TIMER_STOP(t_vec_aypx);    PetscFunctionReturn(0);} // VecAYPX

Pointer<SAMRAIVectorReal<NDIM, double> > PETScNewtonKrylovSolver::getFunctionVector() const{    Vec petsc_f;    int ierr = SNESGetFunction(d_petsc_snes, &petsc_f, NULL, NULL);    IBTK_CHKERRQ(ierr);    return PETScSAMRAIVectorReal::getSAMRAIVector(petsc_f);} // getFunctionVector

Pointer<SAMRAIVectorReal<NDIM, double> > PETScSNESJacobianJOWrapper::getBaseVector() const{    Vec petsc_x;    int ierr = SNESGetSolution(d_petsc_snes, &petsc_x);    IBTK_CHKERRQ(ierr);    return PETScSAMRAIVectorReal::getSAMRAIVector(petsc_x);} // getBaseVector

PetscErrorCode VecSetRandom_SAMRAI(Vec x, PetscRandom rctx){    IBTK_TIMER_START(t_vec_set_random);#if !defined(NDEBUG)    TBOX_ASSERT(x);#endif    PetscScalar lo, hi;    int ierr;    ierr = PetscRandomGetInterval(rctx, &lo, &hi);    IBTK_CHKERRQ(ierr);    PSVR_CAST2(x)->setRandomValues(hi - lo, lo);    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(x));    IBTK_CHKERRQ(ierr);    IBTK_TIMER_STOP(t_vec_set_random);    PetscFunctionReturn(0);} // VecSetRandom

PetscErrorCodeVecMAXPY_SAMRAI(    Vec y,    PetscInt nv,    const PetscScalar* alpha,    Vec* x){    IBTK_TIMER_START(t_vec_maxpy);#ifdef DEBUG_CHECK_ASSERTIONS    TBOX_ASSERT(y != PETSC_NULL);    for (PetscInt i = 0; i < nv; ++i)    {        TBOX_ASSERT(x[i] != PETSC_NULL);    }#endif    static const bool interior_only = false;    for (PetscInt i = 0; i < nv; ++i)    {        if (MathUtilities<double>::equalEps(alpha[i],1.0))        {            PSVR_CAST2(y)->add(PSVR_CAST2(x[i]), PSVR_CAST2(y), interior_only);        }        else if (MathUtilities<double>::equalEps(alpha[i],-1.0))        {            PSVR_CAST2(y)->subtract(PSVR_CAST2(y), PSVR_CAST2(x[i]), interior_only);        }        else        {            PSVR_CAST2(y)->axpy(alpha[i], PSVR_CAST2(x[i]), PSVR_CAST2(y), interior_only);        }    }    int ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(y)); IBTK_CHKERRQ(ierr);    IBTK_TIMER_STOP(t_vec_maxpy);    PetscFunctionReturn(0);}// VecMAXPY

PetscErrorCodeVecWAXPY_SAMRAI(    Vec w,    PetscScalar alpha,    Vec x,    Vec y){    IBTK_TIMER_START(t_vec_waxpy);#ifdef DEBUG_CHECK_ASSERTIONS    TBOX_ASSERT(x != PETSC_NULL);    TBOX_ASSERT(y != PETSC_NULL);    TBOX_ASSERT(w != PETSC_NULL);#endif    static const bool interior_only = false;    if (MathUtilities<double>::equalEps(alpha,1.0))    {        PSVR_CAST2(w)->add(PSVR_CAST2(x), PSVR_CAST2(y), interior_only);    }    else if (MathUtilities<double>::equalEps(alpha,-1.0))    {        PSVR_CAST2(w)->subtract(PSVR_CAST2(y), PSVR_CAST2(x), interior_only);    }    else    {        PSVR_CAST2(w)->axpy(alpha, PSVR_CAST2(x), PSVR_CAST2(y), interior_only);    }    int ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(w)); IBTK_CHKERRQ(ierr);    IBTK_TIMER_STOP(t_vec_waxpy);    PetscFunctionReturn(0);}// VecWAXPY

PetscErrorCode VecSwap_SAMRAI(Vec x, Vec y){    IBTK_TIMER_START(t_vec_swap);#if !defined(NDEBUG)    TBOX_ASSERT(x);    TBOX_ASSERT(y);#endif    PSVR_CAST2(x)->swapVectors(PSVR_CAST2(y));    int ierr;    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(x));    IBTK_CHKERRQ(ierr);    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(y));    IBTK_CHKERRQ(ierr);    IBTK_TIMER_STOP(t_vec_swap);    PetscFunctionReturn(0);} // VecSwap

LData::LData(const std::string& name,             const unsigned int num_local_nodes,             const unsigned int depth,             const std::vector<int>& nonlocal_petsc_indices)    : d_name(name), d_global_node_count(0), d_local_node_count(0), d_ghost_node_count(0), d_depth(depth),      d_nonlocal_petsc_indices(nonlocal_petsc_indices), d_global_vec(NULL), d_managing_petsc_vec(true), d_array(NULL),      d_boost_array(NULL), d_boost_local_array(NULL), d_boost_vec_array(NULL), d_boost_local_vec_array(NULL),      d_ghosted_local_vec(NULL), d_ghosted_local_array(NULL), d_boost_ghosted_local_array(NULL),      d_boost_vec_ghosted_local_array(NULL){    // Create the PETSc Vec that provides storage for the Lagrangian data.    int ierr;    if (d_depth == 1)    {        ierr = VecCreateGhost(PETSC_COMM_WORLD,                              num_local_nodes,                              PETSC_DECIDE,                              static_cast<int>(d_nonlocal_petsc_indices.size()),                              d_nonlocal_petsc_indices.empty() ? NULL : &d_nonlocal_petsc_indices[0],                              &d_global_vec);        IBTK_CHKERRQ(ierr);    }    else    {        ierr = VecCreateGhostBlock(PETSC_COMM_WORLD,                                   d_depth,                                   d_depth * num_local_nodes,                                   PETSC_DECIDE,                                   static_cast<int>(d_nonlocal_petsc_indices.size()),                                   d_nonlocal_petsc_indices.empty() ? NULL : &d_nonlocal_petsc_indices[0],                                   &d_global_vec);        IBTK_CHKERRQ(ierr);    }    int global_node_count;    ierr = VecGetSize(d_global_vec, &global_node_count);    IBTK_CHKERRQ(ierr);#if !defined(NDEBUG)    TBOX_ASSERT(global_node_count >= 0);#endif    d_global_node_count = global_node_count;    d_global_node_count /= d_depth;    d_local_node_count = num_local_nodes;    d_ghost_node_count = static_cast<int>(d_nonlocal_petsc_indices.size());    return;} // LData

void GeneralizedIBMethod::computeLagrangianForce(const double data_time){    IBMethod::computeLagrangianForce(data_time);    const int coarsest_ln = 0;    const int finest_ln = d_hierarchy->getFinestLevelNumber();    int ierr;    std::vector<Pointer<LData> >* F_data = NULL;    std::vector<Pointer<LData> >* N_data = NULL;    std::vector<Pointer<LData> >* X_data = NULL;    std::vector<Pointer<LData> >* D_data = NULL;    if (MathUtilities<double>::equalEps(data_time, d_current_time))    {        d_F_current_needs_ghost_fill = true;        d_N_current_needs_ghost_fill = true;        F_data = &d_F_current_data;        N_data = &d_N_current_data;        X_data = &d_X_current_data;        D_data = &d_D_current_data;    }    else if (MathUtilities<double>::equalEps(data_time, d_half_time))    {        TBOX_ERROR(d_object_name << "::computeLagrangianForce():/n"                                 << "  time-stepping type MIDPOINT_RULE not supported by "                                    "class GeneralizedIBMethod;/n"                                 << "  use TRAPEZOIDAL_RULE instead./n");    }    else if (MathUtilities<double>::equalEps(data_time, d_new_time))    {        d_F_new_needs_ghost_fill = true;        d_N_new_needs_ghost_fill = true;        F_data = &d_F_new_data;        N_data = &d_N_new_data;        X_data = &d_X_new_data;        D_data = &d_D_new_data;    }    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)    {        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;        ierr = VecSet((*N_data)[ln]->getVec(), 0.0);        IBTK_CHKERRQ(ierr);        if (d_ib_force_and_torque_fcn)        {            d_ib_force_and_torque_fcn->computeLagrangianForceAndTorque((*F_data)[ln],                                                                       (*N_data)[ln],                                                                       (*X_data)[ln],                                                                       (*D_data)[ln],                                                                       d_hierarchy,                                                                       ln,                                                                       data_time,                                                                       d_l_data_manager);        }    }    resetAnchorPointValues(*F_data, coarsest_ln, finest_ln);    resetAnchorPointValues(*N_data, coarsest_ln, finest_ln);    return;} // computeLagrangianForce

PetscErrorCode IBImplicitStaggeredHierarchyIntegrator::compositeIBJacobianSetup(    SNES /*snes*/,    Vec x,    Mat* A,    Mat* /*B*/,    MatStructure* /*mat_structure*/){    PetscErrorCode ierr;    ierr = MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY);    IBTK_CHKERRQ(ierr);    ierr = MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY);    IBTK_CHKERRQ(ierr);    Vec* component_sol_vecs;    ierr = VecMultiVecGetSubVecs(x, &component_sol_vecs);    IBTK_CHKERRQ(ierr);    Vec X = component_sol_vecs[1];    d_ib_implicit_ops->setLinearizedPosition(X);    return 0;} // compositeIBJacobianSetup

LData::~LData(){    restoreArrays();    if (d_managing_petsc_vec)    {        const int ierr = VecDestroy(&d_global_vec);        IBTK_CHKERRQ(ierr);    }    return;} // ~LData

voidVCSCViscousPETScLevelSolver::initializeSolverStateSpecialized(const SAMRAIVectorReal<NDIM, double>& x,                                                              const SAMRAIVectorReal<NDIM, double>& /*b*/){    // Allocate DOF index data.    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();    const int x_idx = x.getComponentDescriptorIndex(0);    Pointer<SideDataFactory<NDIM, double> > x_fac = var_db->getPatchDescriptor()->getPatchDataFactory(x_idx);    const int depth = x_fac->getDefaultDepth();    Pointer<SideDataFactory<NDIM, int> > dof_index_fac =        var_db->getPatchDescriptor()->getPatchDataFactory(d_dof_index_idx);    dof_index_fac->setDefaultDepth(depth);    if (!d_level->checkAllocated(d_dof_index_idx)) d_level->allocatePatchData(d_dof_index_idx);    PETScVecUtilities::constructPatchLevelDOFIndices(d_num_dofs_per_proc, d_dof_index_idx, d_level);    // Setup PETSc objects.    int ierr;    const int mpi_rank = SAMRAI_MPI::getRank();    ierr = VecCreateMPI(PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_x);    IBTK_CHKERRQ(ierr);    ierr = VecCreateMPI(PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_b);    IBTK_CHKERRQ(ierr);    const double alpha = 1.0;    const double beta = 1.0;    PETScMatUtilities::constructPatchLevelVCSCViscousOp(d_petsc_mat,                                                        d_poisson_spec,                                                        alpha,                                                        beta,                                                        d_bc_coefs,                                                        d_solution_time,                                                        d_num_dofs_per_proc,                                                        d_dof_index_idx,                                                        d_level,                                                        d_mu_interp_type);    d_petsc_pc = d_petsc_mat;    // Setup SAMRAI communication objects.    d_data_synch_sched = PETScVecUtilities::constructDataSynchSchedule(x_idx, d_level);    d_ghost_fill_sched = PETScVecUtilities::constructGhostFillSchedule(x_idx, d_level);    return;} // initializeSolverStateSpecialized

PetscErrorCodeIBImplicitStaggeredHierarchyIntegrator::compositeIBJacobianApply_SAMRAI(Mat A, Vec x, Vec f){    PetscErrorCode ierr;    void* ctx;    ierr = MatShellGetContext(A, &ctx);    IBTK_CHKERRQ(ierr);    IBImplicitStaggeredHierarchyIntegrator* ib_integrator =        static_cast<IBImplicitStaggeredHierarchyIntegrator*>(ctx);    return ib_integrator->compositeIBJacobianApply(x, f);} // compositeIBJacobianApply_SAMRAI

void PETScSNESFunctionGOWrapper::apply(SAMRAIVectorReal<NDIM, double>& x, SAMRAIVectorReal<NDIM, double>& y){    if (!d_is_initialized) initializeOperatorState(x, y);    // Update the PETSc Vec wrappers.    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_x, Pointer<SAMRAIVectorReal<NDIM, double> >(&x, false));    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_y, Pointer<SAMRAIVectorReal<NDIM, double> >(&y, false));    // Apply the operator.    int ierr = d_petsc_snes_form_func(d_petsc_snes, d_petsc_x, d_petsc_y, d_petsc_snes_func_ctx);    IBTK_CHKERRQ(ierr);    return;} // apply

PetscErrorCode VecSet_SAMRAI(Vec x, PetscScalar alpha){    IBTK_TIMER_START(t_vec_set);#if !defined(NDEBUG)    TBOX_ASSERT(x);#endif    static const bool interior_only = false;    PSVR_CAST2(x)->setToScalar(alpha, interior_only);    int ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(x));    IBTK_CHKERRQ(ierr);    IBTK_TIMER_STOP(t_vec_set);    PetscFunctionReturn(0);} // VecSet

void PETScSNESFunctionGOWrapper::initializeOperatorState(const SAMRAIVectorReal<NDIM, double>& in,                                                         const SAMRAIVectorReal<NDIM, double>& out){    if (d_is_initialized) deallocateOperatorState();    d_x = in.cloneVector("");    d_y = out.cloneVector("");    MPI_Comm comm;    int ierr = PetscObjectGetComm(reinterpret_cast<PetscObject>(d_petsc_snes), &comm);    IBTK_CHKERRQ(ierr);    d_petsc_x = PETScSAMRAIVectorReal::createPETScVector(d_x, comm);    d_petsc_y = PETScSAMRAIVectorReal::createPETScVector(d_y, comm);    d_is_initialized = true;    return;} // initializeOperatorState

boolPETScPCLSWrapper::solveSystem(SAMRAIVectorReal<NDIM, double>& x, SAMRAIVectorReal<NDIM, double>& b){    if (!d_is_initialized) initializeSolverState(x, b);    // Update the PETSc Vec wrappers.    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_x, Pointer<SAMRAIVectorReal<NDIM, double> >(&x, false));    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_b, Pointer<SAMRAIVectorReal<NDIM, double> >(&b, false));    // Apply the preconditioner.    int ierr = PCApply(d_petsc_pc, d_petsc_x, d_petsc_b);    IBTK_CHKERRQ(ierr);    return true;} // solveSystem

PetscErrorCode VecCopy_SAMRAI(Vec x, Vec y){    IBTK_TIMER_START(t_vec_copy);#if !defined(NDEBUG)    TBOX_ASSERT(x);    TBOX_ASSERT(y);#endif    static const bool interior_only = false;    PSVR_CAST2(y)->copyVector(PSVR_CAST2(x), interior_only);    int ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(y));    IBTK_CHKERRQ(ierr);    IBTK_TIMER_STOP(t_vec_copy);    PetscFunctionReturn(0);} // VecCopy