C++ Patch::getBox方法代码示例

BoundaryBox<NDIM>
PhysicalBoundaryUtilities::trimBoundaryCodim1Box(const BoundaryBox<NDIM>& bdry_box,
                                                 const Patch<NDIM>& patch)
{
#if !defined(NDEBUG)
    TBOX_ASSERT(bdry_box.getBoundaryType() == 1);
#endif
    // Trim a boundary box so it does not stick out past the corners of a patch.
    const Box<NDIM>& b_box = bdry_box.getBox();
    const Box<NDIM>& patch_box = patch.getBox();
    const unsigned int bdry_normal_axis = bdry_box.getLocationIndex() / 2;

    Box<NDIM> trimmed_b_box = b_box;
    for (unsigned int d = 0; d < NDIM; ++d)
    {
        if (d != bdry_normal_axis)
        {
            trimmed_b_box.lower()[d] = std::max(b_box.lower()[d], patch_box.lower()[d]);
            trimmed_b_box.upper()[d] = std::min(b_box.upper()[d], patch_box.upper()[d]);
        }
    }
    const BoundaryBox<NDIM> trimmed_bdry_box(
        trimmed_b_box, bdry_box.getBoundaryType(), bdry_box.getLocationIndex());
    return trimmed_bdry_box;
} // trimBoundaryCodim1Box

void CartCellRobinPhysBdryOp::fillGhostCellValuesCodim2(
    const int patch_data_idx,
    const Array<BoundaryBox<NDIM> >& physical_codim2_boxes,
    const IntVector<NDIM>& ghost_width_to_fill,
    const Patch<NDIM>& patch,
    const bool adjoint_op)
{
    const int n_physical_codim2_boxes = physical_codim2_boxes.size();
    if (n_physical_codim2_boxes == 0) return;

    const Box<NDIM>& patch_box = patch.getBox();
    Pointer<CartesianPatchGeometry<NDIM> > pgeom = patch.getPatchGeometry();
    Pointer<CellData<NDIM, double> > patch_data = patch.getPatchData(patch_data_idx);
    const int patch_data_depth = patch_data->getDepth();
    const int patch_data_gcw = (patch_data->getGhostCellWidth()).max();
#if !defined(NDEBUG)
    if (patch_data_gcw != (patch_data->getGhostCellWidth()).min())
    {
        TBOX_ERROR(
            "CartCellRobinPhysBdryOp::fillGhostCellValuesCodim2():\n"
            "  patch data for patch data index "
            << patch_data_idx << " does not have uniform ghost cell widths." << std::endl);
    }
#endif
    const IntVector<NDIM> gcw_to_fill =
        IntVector<NDIM>::min(patch_data->getGhostCellWidth(), ghost_width_to_fill);

    for (int n = 0; n < n_physical_codim2_boxes; ++n)
    {
        const BoundaryBox<NDIM>& bdry_box = physical_codim2_boxes[n];
        const unsigned int location_index = bdry_box.getLocationIndex();
        const Box<NDIM> bc_fill_box =
            pgeom->getBoundaryFillBox(bdry_box, patch_box, gcw_to_fill);
        for (int d = 0; d < patch_data_depth; ++d)
        {
            CC_ROBIN_PHYS_BDRY_OP_2_FC(patch_data->getPointer(d),
                                       patch_data_gcw,
                                       location_index,
                                       patch_box.lower(0),
                                       patch_box.upper(0),
                                       patch_box.lower(1),
                                       patch_box.upper(1),
#if (NDIM == 3)
                                       patch_box.lower(2),
                                       patch_box.upper(2),
#endif
                                       bc_fill_box.lower(0),
                                       bc_fill_box.upper(0),
                                       bc_fill_box.lower(1),
                                       bc_fill_box.upper(1),
#if (NDIM == 3)
                                       bc_fill_box.lower(2),
                                       bc_fill_box.upper(2),
#endif
                                       adjoint_op ? 1 : 0);
        }
    }
    return;
} // fillGhostCellValuesCodim2

void
VecCellRefineAdapter::refine(
    Patch<NDIM>& fine,
    const Patch<NDIM>& coarse,
    const int dst_component,
    const int src_component,
    const Box<NDIM>& fine_box,
    const IntVector<NDIM>& ratio) const
{
    Pointer<VecCellData<double> > dst_data = fine  .getPatchData(dst_component);
    Pointer<VecCellData<double> > src_data = coarse.getPatchData(src_component);
#ifdef DEBUG_CHECK_ASSERTIONS
    TBOX_ASSERT(!dst_data.isNull());
    TBOX_ASSERT(!src_data.isNull());
#endif
    // Create "dummy" patches.
    Patch<NDIM> fine_cell(fine.getBox(), fine.getPatchDescriptor());
    fine_cell.setPatchGeometry(fine.getPatchGeometry());
    fine_cell.setPatchInHierarchy(fine.inHierarchy());
    fine_cell.setPatchLevelNumber(fine.getPatchLevelNumber());
    fine_cell.setPatchNumber(fine.getPatchNumber());

    Patch<NDIM> coarse_cell(coarse.getBox(), coarse.getPatchDescriptor());
    coarse_cell.setPatchGeometry(coarse.getPatchGeometry());
    coarse_cell.setPatchInHierarchy(coarse.inHierarchy());
    coarse_cell.setPatchLevelNumber(coarse.getPatchLevelNumber());
    coarse_cell.setPatchNumber(coarse.getPatchNumber());

    // Make copies of the dst and src data.
    CellData<NDIM,double> dst_cell_data(dst_data->getBox(), dst_data->getDepth(), dst_data->getGhostCellWidth());
    dst_data->copy2(dst_cell_data);
    fine_cell.allocatePatchData(dst_component);
    fine_cell.setPatchData(dst_component,Pointer<PatchData<NDIM> >(&dst_cell_data,false));

    CellData<NDIM,double> src_cell_data(src_data->getBox(), src_data->getDepth(), src_data->getGhostCellWidth());
    src_data->copy2(src_cell_data);
    coarse_cell.allocatePatchData(src_component);
    coarse_cell.setPatchData(src_component,Pointer<PatchData<NDIM> >(&src_cell_data,false));

    // Refine data from the coarse dummy patch to the fine dummy patch.
    d_cell_refine_op->refine(fine_cell, coarse_cell, dst_component, src_component, fine_box, ratio);

    // Copy the result into the VecCellData.
    dst_data->copy(dst_cell_data);
    return;
}// refine

void LMarkerRefine::refine(Patch<NDIM>& fine,
                           const Patch<NDIM>& coarse,
                           const int dst_component,
                           const int src_component,
                           const Box<NDIM>& fine_box,
                           const IntVector<NDIM>& ratio) const
{
    Pointer<LMarkerSetData> dst_mark_data = fine.getPatchData(dst_component);
    Pointer<LMarkerSetData> src_mark_data = coarse.getPatchData(src_component);

    const Box<NDIM>& fine_patch_box = fine.getBox();
    const Pointer<CartesianPatchGeometry<NDIM> > fine_patch_geom = fine.getPatchGeometry();
    const Index<NDIM>& fine_patch_lower = fine_patch_box.lower();
    const Index<NDIM>& fine_patch_upper = fine_patch_box.upper();
    const double* const fine_patchXLower = fine_patch_geom->getXLower();
    const double* const fine_patchXUpper = fine_patch_geom->getXUpper();
    const double* const fine_patchDx = fine_patch_geom->getDx();

    const Pointer<CartesianPatchGeometry<NDIM> > coarse_patch_geom = coarse.getPatchGeometry();
    const double* const coarse_patchDx = coarse_patch_geom->getDx();

    const Box<NDIM> coarse_box = Box<NDIM>::coarsen(fine_box, ratio);
    for (LMarkerSetData::SetIterator it(*src_mark_data); it; it++)
    {
        const Index<NDIM>& coarse_i = it.getIndex();
        if (coarse_box.contains(coarse_i))
        {
            const LMarkerSet& coarse_mark_set = it();
            for (LMarkerSet::const_iterator cit = coarse_mark_set.begin(); cit != coarse_mark_set.end(); ++cit)
            {
                const LMarkerSet::value_type& coarse_mark = *cit;
                const Point& X = coarse_mark->getPosition();
                const IntVector<NDIM>& offset = coarse_mark->getPeriodicOffset();
                boost::array<double, NDIM> X_shifted;
                for (unsigned int d = 0; d < NDIM; ++d)
                {
                    X_shifted[d] = X[d] + static_cast<double>(offset(d)) * coarse_patchDx[d];
                }

                const Index<NDIM> fine_i = IndexUtilities::getCellIndex(
                    X_shifted, fine_patchXLower, fine_patchXUpper, fine_patchDx, fine_patch_lower, fine_patch_upper);
                if (fine_box.contains(fine_i))
                {
                    if (!dst_mark_data->isElement(fine_i))
                    {
                        dst_mark_data->appendItemPointer(fine_i, new LMarkerSet());
                    }
                    LMarkerSet& fine_mark_set = *(dst_mark_data->getItem(fine_i));
                    fine_mark_set.push_back(coarse_mark);
                }
            }
        }
    }
    return;
} // refine

void CartCellDoubleLinearCFInterpolation::computeNormalExtension(
    Patch<NDIM>& patch,
    const IntVector<NDIM>& ratio,
    const IntVector<NDIM>& /*ghost_width_to_fill*/)
{
#if !defined(NDEBUG)
    TBOX_ASSERT(d_hierarchy);
    TBOX_ASSERT(!d_consistent_type_2_bdry);
    TBOX_ASSERT(ratio.min() == ratio.max());
#endif
    // Ensure that the fine patch is located on the expected destination level;
    // if not, we are not guaranteed to have appropriate coarse-fine interface
    // boundary box information.
    if (!patch.inHierarchy()) return;

    // Get the co-dimension 1 cf boundary boxes.
    const int patch_num = patch.getPatchNumber();
    const int patch_level_num = patch.getPatchLevelNumber();
#if !defined(NDEBUG)
    Pointer<PatchLevel<NDIM> > level = d_hierarchy->getPatchLevel(patch_level_num);
    TBOX_ASSERT(&patch == level->getPatch(patch_num).getPointer());
#endif
    const Array<BoundaryBox<NDIM> >& cf_bdry_codim1_boxes =
        d_cf_boundary[patch_level_num]->getBoundaries(patch_num, 1);
    const int n_cf_bdry_codim1_boxes = cf_bdry_codim1_boxes.size();

    // Check to see if there are any co-dimension 1 coarse-fine boundary boxes
    // associated with the patch; if not, there is nothing to do.
    if (n_cf_bdry_codim1_boxes == 0) return;

    // Get the patch data.
    for (std::set<int>::const_iterator cit = d_patch_data_indices.begin();
            cit != d_patch_data_indices.end();
            ++cit)
    {
        const int& patch_data_index = *cit;
        Pointer<CellData<NDIM, double> > data = patch.getPatchData(patch_data_index);
#if !defined(NDEBUG)
        TBOX_ASSERT(data);
#endif
        const int U_ghosts = (data->getGhostCellWidth()).max();
#if !defined(NDEBUG)
        if (U_ghosts != (data->getGhostCellWidth()).min())
        {
            TBOX_ERROR("CartCellDoubleLinearCFInterpolation::computeNormalExtension():\n"
                       << "   patch data does not have uniform ghost cell widths"
                       << std::endl);
        }
#endif
        const int data_depth = data->getDepth();
        const IntVector<NDIM> ghost_width_to_fill = GHOST_WIDTH_TO_FILL;
        Pointer<CartesianPatchGeometry<NDIM> > pgeom = patch.getPatchGeometry();
        const Box<NDIM>& patch_box = patch.getBox();
        for (int k = 0; k < n_cf_bdry_codim1_boxes; ++k)
        {
            const BoundaryBox<NDIM>& bdry_box = cf_bdry_codim1_boxes[k];
            const Box<NDIM> bc_fill_box =
                pgeom->getBoundaryFillBox(bdry_box, patch_box, ghost_width_to_fill);
            const unsigned int location_index = bdry_box.getLocationIndex();
            for (int depth = 0; depth < data_depth; ++depth)
            {
                double* const U = data->getPointer(depth);
                const int r = ratio.min();
                CC_LINEAR_NORMAL_INTERPOLATION_FC(U,
                                                  U_ghosts,
                                                  patch_box.lower(0),
                                                  patch_box.upper(0),
                                                  patch_box.lower(1),
                                                  patch_box.upper(1),
#if (NDIM == 3)
                                                  patch_box.lower(2),
                                                  patch_box.upper(2),
#endif
                                                  location_index,
                                                  r,
                                                  bc_fill_box.lower(),
                                                  bc_fill_box.upper());
            }
        }
    }
    return;
} // computeNormalExtension

void
CartSideDoubleDivPreservingRefine::postprocessRefine(Patch<NDIM>& fine,
                                                     const Patch<NDIM>& coarse,
                                                     const Box<NDIM>& unrestricted_fine_box,
                                                     const IntVector<NDIM>& ratio)
{
    // NOTE: This operator cannot fill the full ghost cell width of the
    // destination data.  We instead restrict the size of the fine box to ensure
    // that we have adequate data to apply the divergence- and curl-preserving
    // corrections.
    const Box<NDIM> fine_box = unrestricted_fine_box * Box<NDIM>::grow(fine.getBox(), 2);

#if !defined(NDEBUG)
    for (int d = 0; d < NDIM; ++d)
    {
        if (ratio(d) % 2 != 0)
        {
            TBOX_ERROR("CartSideDoubleDivPreservingRefine::postprocessRefine():\n"
                       << "  refinement ratio must be a power of 2 for divergence- and "
                          "curl-preserving refinement operator."
                       << std::endl);
        }
    }
#endif

    Pointer<SideData<NDIM, double> > fdata = fine.getPatchData(d_u_dst_idx);
#if !defined(NDEBUG)
    TBOX_ASSERT(fdata);
#endif
    const int fdata_ghosts = fdata->getGhostCellWidth().max();
#if !defined(NDEBUG)
    TBOX_ASSERT(fdata_ghosts == fdata->getGhostCellWidth().min());
#endif
    const int fdata_depth = fdata->getDepth();

    Pointer<SideData<NDIM, double> > cdata = coarse.getPatchData(d_u_dst_idx);
#if !defined(NDEBUG)
    TBOX_ASSERT(cdata);
    const int cdata_ghosts = cdata->getGhostCellWidth().max();
    TBOX_ASSERT(cdata_ghosts == cdata->getGhostCellWidth().min());
    const int cdata_depth = cdata->getDepth();
    TBOX_ASSERT(cdata_depth == fdata_depth);
#endif

    if (ratio == IntVector<NDIM>(2))
    {
        // Perform (limited) conservative prolongation of the coarse grid data.
        d_refine_op->refine(fine, coarse, d_u_dst_idx, d_u_dst_idx, fine_box, ratio);

        Pointer<SideData<NDIM, double> > u_src_data = fine.getPatchData(d_u_src_idx);
        Pointer<SideData<NDIM, double> > indicator_data = fine.getPatchData(d_indicator_idx);

        // Ensure that we do not modify any of the data from the old level by
        // setting the value of the fine grid data to equal u_src wherever the
        // indicator data equals "1".
        if (u_src_data && indicator_data)
        {
            for (unsigned int axis = 0; axis < NDIM; ++axis)
            {
                for (Box<NDIM>::Iterator b(SideGeometry<NDIM>::toSideBox(fine_box, axis)); b; b++)
                {
                    const Index<NDIM>& i = b();
                    const SideIndex<NDIM> i_s(i, axis, 0);
                    if (std::abs((*indicator_data)(i_s)-1.0) < 1.0e-12)
                    {
                        for (int depth = 0; depth < fdata_depth; ++depth)
                        {
                            (*fdata)(i_s, depth) = (*u_src_data)(i_s, depth);
                        }
                    }
                }
            }
        }

        // Reinterpolate data in the normal direction in the newly refined part
        // of the level wherever the indicator data does NOT equal "1".  Notice
        // that this loop actually modifies only data that is NOT covered by an
        // overlying coarse grid cell face.
        if (indicator_data)
        {
            for (unsigned int axis = 0; axis < NDIM; ++axis)
            {
                for (Box<NDIM>::Iterator b(SideGeometry<NDIM>::toSideBox(fine_box, axis)); b; b++)
                {
                    const Index<NDIM>& i = b();
                    const SideIndex<NDIM> i_s(i, axis, 0);
                    if (!(std::abs((*indicator_data)(i_s)-1.0) < 1.0e-12))
                    {
                        const Index<NDIM> i_coarse_lower = IndexUtilities::coarsen(i, ratio);
                        const Index<NDIM> i_lower = IndexUtilities::refine(i_coarse_lower, ratio);
                        if (i(axis) == i_lower(axis)) continue;

                        Index<NDIM> i_coarse_upper = i_coarse_lower;
                        i_coarse_upper(axis) += 1;
                        const Index<NDIM> i_upper = IndexUtilities::refine(i_coarse_upper, ratio);

                        const double w1 =
                            static_cast<double>(i(axis) - i_lower(axis)) / static_cast<double>(ratio(axis));
                        const double w0 = 1.0 - w1;

//.........这里部分代码省略.........

void
CartCellDoubleBoundsPreservingConservativeLinearRefine::refine(Patch<NDIM>& fine,
                                                               const Patch<NDIM>& coarse,
                                                               const int dst_component,
                                                               const int src_component,
                                                               const Box<NDIM>& fine_box,
                                                               const IntVector<NDIM>& ratio)
    const
{
    // Determine the box over which we can apply the bounds-preserving
    // correction, and construct a list of boxes that will not be corrected.
    bool empty_correction_box = false;
    Box<NDIM> correction_box = Box<NDIM>::refine(Box<NDIM>::coarsen(fine_box, ratio), ratio);
    for (unsigned int axis = 0; axis < NDIM; ++axis)
    {
        int& lower = correction_box.lower()(axis);
        while (lower < fine_box.lower()(axis))
        {
            lower += ratio(axis);
        }

        int& upper = correction_box.upper()(axis);
        while (upper > fine_box.upper()(axis))
        {
            upper -= ratio(axis);
        }

        if (lower >= upper)
        {
            empty_correction_box = true;
        }
    }
    const Box<NDIM> coarse_correction_box = Box<NDIM>::coarsen(correction_box, ratio);

    BoxList<NDIM> uncorrected_boxes(fine_box);
    if (!empty_correction_box)
    {
        uncorrected_boxes.removeIntersections(correction_box);
    }

    // Employ limited conservative interpolation to prolong data on the
    // correction box.
    d_conservative_linear_refine_op.refine(
        fine, coarse, dst_component, src_component, correction_box, ratio);

    // Employ constant interpolation to prolong data on the rest of the fine
    // box.
    for (BoxList<NDIM>::Iterator b(uncorrected_boxes); b; b++)
    {
        d_constant_refine_op.refine(fine, coarse, dst_component, src_component, b(), ratio);
    }

    // There is nothing left to do if the correction box is empty.
    if (empty_correction_box) return;

    // Correct the data within the correction box.
    Pointer<CellData<NDIM, double> > fdata = fine.getPatchData(dst_component);
    Pointer<CellData<NDIM, double> > cdata = coarse.getPatchData(src_component);
#if !defined(NDEBUG)
    TBOX_ASSERT(fdata);
    TBOX_ASSERT(cdata);
    TBOX_ASSERT(fdata->getDepth() == cdata->getDepth());
#endif
    const int data_depth = fdata->getDepth();
    const Box<NDIM>& patch_box_crse = coarse.getBox();
    const Index<NDIM>& patch_lower_crse = patch_box_crse.lower();
    const Index<NDIM>& patch_upper_crse = patch_box_crse.upper();
    Pointer<CartesianPatchGeometry<NDIM> > pgeom_crse = coarse.getPatchGeometry();
    for (int depth = 0; depth < data_depth; ++depth)
    {
        for (Box<NDIM>::Iterator b(coarse_correction_box); b; b++)
        {
            const Index<NDIM>& i_crse = b();
            const Index<NDIM> i_fine = i_crse * ratio;

            // Determine the lower/upper bounds.
            Box<NDIM> stencil_box_crse(i_crse, i_crse);
            for (unsigned int axis = 0; axis < NDIM; ++axis)
            {
                if (i_crse(axis) > patch_lower_crse(axis) ||
                    !pgeom_crse->getTouchesRegularBoundary(axis, 0))
                {
                    stencil_box_crse.growLower(axis, 1);
                }
                if (i_crse(axis) < patch_upper_crse(axis) ||
                    !pgeom_crse->getTouchesRegularBoundary(axis, 1))
                {
                    stencil_box_crse.growUpper(axis, 1);
                }
            }

            double l = std::numeric_limits<double>::max();
            double u = -(l - std::numeric_limits<double>::epsilon());
            for (Box<NDIM>::Iterator b(stencil_box_crse); b; b++)
            {
                const double& m = (*cdata)(b(), depth);
                l = std::min(l, m);
                u = std::max(u, m);
            }

//.........这里部分代码省略.........

void CartCellRobinPhysBdryOp::fillGhostCellValuesCodim1(
    const int patch_data_idx,
    const Array<BoundaryBox<NDIM> >& physical_codim1_boxes,
    const double fill_time,
    const IntVector<NDIM>& ghost_width_to_fill,
    Patch<NDIM>& patch,
    const bool adjoint_op)
{
    const int n_physical_codim1_boxes = physical_codim1_boxes.size();
    if (n_physical_codim1_boxes == 0) return;

    const Box<NDIM>& patch_box = patch.getBox();
    Pointer<CartesianPatchGeometry<NDIM> > pgeom = patch.getPatchGeometry();
    const double* const dx = pgeom->getDx();
    Pointer<CellData<NDIM, double> > patch_data = patch.getPatchData(patch_data_idx);
    const int patch_data_depth = patch_data->getDepth();
    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
    Pointer<Variable<NDIM> > var;
    var_db->mapIndexToVariable(patch_data_idx, var);
    const int patch_data_gcw = (patch_data->getGhostCellWidth()).max();
#if !defined(NDEBUG)
    if (patch_data_gcw != (patch_data->getGhostCellWidth()).min())
    {
        TBOX_ERROR(
            "CartCellRobinPhysBdryOp::fillGhostCellValuesCodim1():\n"
            "  patch data for patch data index "
            << patch_data_idx << " does not have uniform ghost cell widths." << std::endl);
    }
#endif
    const IntVector<NDIM> gcw_to_fill =
        IntVector<NDIM>::min(patch_data->getGhostCellWidth(), ghost_width_to_fill);

    // Set the boundary condition coefficients and then set the ghost cell
    // values.
    for (int n = 0; n < n_physical_codim1_boxes; ++n)
    {
        const BoundaryBox<NDIM>& bdry_box = physical_codim1_boxes[n];
        const unsigned int location_index = bdry_box.getLocationIndex();
        const Box<NDIM> bc_fill_box =
            pgeom->getBoundaryFillBox(bdry_box, patch_box, gcw_to_fill);
        const BoundaryBox<NDIM> trimmed_bdry_box(bdry_box.getBox() * bc_fill_box,
                                                 bdry_box.getBoundaryType(),
                                                 bdry_box.getLocationIndex());
        const Box<NDIM> bc_coef_box =
            PhysicalBoundaryUtilities::makeSideBoundaryCodim1Box(trimmed_bdry_box);
        Pointer<ArrayData<NDIM, double> > acoef_data =
            new ArrayData<NDIM, double>(bc_coef_box, 1);
        Pointer<ArrayData<NDIM, double> > bcoef_data =
            new ArrayData<NDIM, double>(bc_coef_box, 1);
        Pointer<ArrayData<NDIM, double> > gcoef_data =
            new ArrayData<NDIM, double>(bc_coef_box, 1);
        for (int d = 0; d < patch_data_depth; ++d)
        {
            RobinBcCoefStrategy<NDIM>* bc_coef = d_bc_coefs[d];
            ExtendedRobinBcCoefStrategy* const extended_bc_coef =
                dynamic_cast<ExtendedRobinBcCoefStrategy*>(bc_coef);
            if (extended_bc_coef)
            {
                extended_bc_coef->setTargetPatchDataIndex(patch_data_idx);
                extended_bc_coef->setHomogeneousBc(d_homogeneous_bc);
            }
            bc_coef->setBcCoefs(
                acoef_data, bcoef_data, gcoef_data, var, patch, trimmed_bdry_box, fill_time);
            if (d_homogeneous_bc && !extended_bc_coef) gcoef_data->fillAll(0.0);
            if (extended_bc_coef) extended_bc_coef->clearTargetPatchDataIndex();
            switch (location_index)
            {
            case 0: // lower x
            case 1: // upper x
                CC_ROBIN_PHYS_BDRY_OP_1_X_FC(patch_data->getPointer(d),
                                             patch_data_gcw,
                                             acoef_data->getPointer(),
                                             bcoef_data->getPointer(),
                                             gcoef_data->getPointer(),
                                             location_index,
                                             patch_box.lower(0),
                                             patch_box.upper(0),
                                             patch_box.lower(1),
                                             patch_box.upper(1),
#if (NDIM == 3)
                                             patch_box.lower(2),
                                             patch_box.upper(2),
#endif
                                             bc_fill_box.lower(1),
                                             bc_fill_box.upper(1),
#if (NDIM == 3)
                                             bc_fill_box.lower(2),
                                             bc_fill_box.upper(2),
#endif
                                             dx,
                                             adjoint_op ? 1 : 0);
                break;
            case 2: // lower y
            case 3: // upper y
                CC_ROBIN_PHYS_BDRY_OP_1_Y_FC(patch_data->getPointer(d),
                                             patch_data_gcw,
                                             acoef_data->getPointer(),
                                             bcoef_data->getPointer(),
                                             gcoef_data->getPointer(),
                                             location_index,
//.........这里部分代码省略.........

void muParserRobinBcCoefs::setBcCoefs(Pointer<ArrayData<NDIM, double> >& acoef_data,
                                      Pointer<ArrayData<NDIM, double> >& bcoef_data,
                                      Pointer<ArrayData<NDIM, double> >& gcoef_data,
                                      const Pointer<Variable<NDIM> >& /*variable*/,
                                      const Patch<NDIM>& patch,
                                      const BoundaryBox<NDIM>& bdry_box,
                                      double fill_time) const
{
    const Box<NDIM>& patch_box = patch.getBox();
    const Index<NDIM>& patch_lower = patch_box.lower();
    Pointer<CartesianPatchGeometry<NDIM> > pgeom = patch.getPatchGeometry();

    const double* const x_lower = pgeom->getXLower();
    const double* const dx = pgeom->getDx();

    // Loop over the boundary box and set the coefficients.
    const unsigned int location_index = bdry_box.getLocationIndex();
    const unsigned int bdry_normal_axis = location_index / 2;
    const Box<NDIM>& bc_coef_box =
        (acoef_data ? acoef_data->getBox() : bcoef_data ? bcoef_data->getBox() : gcoef_data ? gcoef_data->getBox() :
                                                                                              Box<NDIM>());
#if !defined(NDEBUG)
    TBOX_ASSERT(!acoef_data || bc_coef_box == acoef_data->getBox());
    TBOX_ASSERT(!bcoef_data || bc_coef_box == bcoef_data->getBox());
    TBOX_ASSERT(!gcoef_data || bc_coef_box == gcoef_data->getBox());
#endif

    const mu::Parser& acoef_parser = d_acoef_parsers[location_index];
    const mu::Parser& bcoef_parser = d_bcoef_parsers[location_index];
    const mu::Parser& gcoef_parser = d_gcoef_parsers[location_index];
    *d_parser_time = fill_time;
    for (Box<NDIM>::Iterator b(bc_coef_box); b; b++)
    {
        const Index<NDIM>& i = b();
        for (unsigned int d = 0; d < NDIM; ++d)
        {
            if (d != bdry_normal_axis)
            {
                (*d_parser_posn)[d] = x_lower[d] + dx[d] * (static_cast<double>(i(d) - patch_lower(d)) + 0.5);
            }
            else
            {
                (*d_parser_posn)[d] = x_lower[d] + dx[d] * (static_cast<double>(i(d) - patch_lower(d)));
            }
        }
        try
        {
            if (acoef_data) (*acoef_data)(i, 0) = acoef_parser.Eval();
            if (bcoef_data) (*bcoef_data)(i, 0) = bcoef_parser.Eval();
            if (gcoef_data) (*gcoef_data)(i, 0) = gcoef_parser.Eval();
        }
        catch (mu::ParserError& e)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  error: " << e.GetMsg() << "\n"
                       << "  in:    " << e.GetExpr() << "\n");
        }
        catch (...)
        {
            TBOX_ERROR("muParserRobinBcCoefs::setDataOnPatch():\n"
                       << "  unrecognized exception generated by muParser library.\n");
        }
    }
    return;
} // setBcCoefs

void
AdvDiffHypPatchOps::conservativeDifferenceOnPatch(
    Patch<NDIM>& patch,
    const double /*time*/,
    const double dt,
    bool /*at_synchronization*/)
{
    const Box<NDIM>& patch_box = patch.getBox();
    const Index<NDIM>& ilower = patch_box.lower();
    const Index<NDIM>& iupper = patch_box.upper();

    const Pointer<CartesianPatchGeometry<NDIM> > patch_geom = patch.getPatchGeometry();
    const double* const dx = patch_geom->getDx();

    for (std::set<Pointer<CellVariable<NDIM,double> > >::const_iterator cit = d_Q_var.begin();
         cit != d_Q_var.end(); ++cit)
    {
        Pointer<CellVariable<NDIM,double> > Q_var = *cit;
        Pointer<FaceVariable<NDIM,double> > u_var = d_Q_u_map[Q_var];

        if (u_var.isNull())
        {
            Pointer<CellData<NDIM,double> > Q_data = patch.getPatchData(Q_var, getDataContext());
            Q_data->fillAll(0.0);
            continue;
        }

        const bool conservation_form = d_Q_difference_form[Q_var] == CONSERVATIVE;
        const bool u_is_div_free = d_u_is_div_free[u_var];

        Pointer<FaceVariable<NDIM,double> > flux_integral_var = d_flux_integral_var[Q_var];
        Pointer<FaceVariable<NDIM,double> > q_integral_var = d_q_integral_var[Q_var];
        Pointer<FaceVariable<NDIM,double> > u_integral_var = d_u_integral_var[u_var];

        Pointer<CellData<NDIM,double> > Q_data = patch.getPatchData(Q_var, getDataContext());
        Pointer<FaceData<NDIM,double> > flux_integral_data =
            (conservation_form
             ? patch.getPatchData(flux_integral_var, getDataContext())
             : Pointer<PatchData<NDIM> >(NULL));
        Pointer<FaceData<NDIM,double> > q_integral_data =
            (!conservation_form || !u_is_div_free
             ? patch.getPatchData(q_integral_var, getDataContext())
             : Pointer<PatchData<NDIM> >(NULL));
        Pointer<FaceData<NDIM,double> > u_integral_data =
            (!conservation_form || !u_is_div_free
             ? patch.getPatchData(u_integral_var, getDataContext())
             : Pointer<PatchData<NDIM> >(NULL));

        const IntVector<NDIM>& Q_data_ghost_cells = Q_data->getGhostCellWidth();
        const IntVector<NDIM>& flux_integral_data_ghost_cells =
            (!flux_integral_data.isNull()
             ? flux_integral_data->getGhostCellWidth()
             : 0);
        const IntVector<NDIM>& q_integral_data_ghost_cells =
            (!q_integral_data.isNull()
             ? q_integral_data->getGhostCellWidth()
             : 0);
        const IntVector<NDIM>& u_integral_data_ghost_cells =
            (!u_integral_data.isNull()
             ? u_integral_data->getGhostCellWidth()
             : 0);

        switch (d_Q_difference_form[Q_var])
        {
            case CONSERVATIVE:
            {
                for (int depth = 0; depth < Q_data->getDepth(); ++depth)
                {
                    if (u_is_div_free)
                    {
#if (NDIM == 2)
                        ADV_DIFF_CONSDIFF_FC(
                            dx,dt,
                            ilower(0),iupper(0),ilower(1),iupper(1),
                            flux_integral_data_ghost_cells(0),flux_integral_data_ghost_cells(1),
                            Q_data_ghost_cells(0),Q_data_ghost_cells(1),
                            flux_integral_data->getPointer(0,depth),
                            flux_integral_data->getPointer(1,depth),
                            Q_data->getPointer(depth));
#endif
#if (NDIM == 3)
                        ADV_DIFF_CONSDIFF_FC(
                            dx,dt,
                            ilower(0),iupper(0),ilower(1),iupper(1),ilower(2),iupper(2),
                            flux_integral_data_ghost_cells(0),flux_integral_data_ghost_cells(1),flux_integral_data_ghost_cells(2),
                            Q_data_ghost_cells(0),Q_data_ghost_cells(1),Q_data_ghost_cells(2),
                            flux_integral_data->getPointer(0,depth),
                            flux_integral_data->getPointer(1,depth),
                            flux_integral_data->getPointer(2,depth),
                            Q_data->getPointer(depth));
#endif
                    }
                    else
                    {
#if (NDIM == 2)
                        ADV_DIFF_CONSDIFFWITHDIVSOURCE_FC(
                            dx,dt,
                            ilower(0),iupper(0),ilower(1),iupper(1),
                            flux_integral_data_ghost_cells(0),flux_integral_data_ghost_cells(1),
                            q_integral_data_ghost_cells(0),q_integral_data_ghost_cells(1),
//.........这里部分代码省略.........

void CartCellDoubleQuadraticRefine::refine(Patch<NDIM>& fine,
                                           const Patch<NDIM>& coarse,
                                           const int dst_component,
                                           const int src_component,
                                           const Box<NDIM>& fine_box,
                                           const IntVector<NDIM>& ratio) const
{
    // Get the patch data.
    Pointer<CellData<NDIM, double> > fdata = fine.getPatchData(dst_component);
    Pointer<CellData<NDIM, double> > cdata = coarse.getPatchData(src_component);
#if !defined(NDEBUG)
    TBOX_ASSERT(fdata);
    TBOX_ASSERT(cdata);
    TBOX_ASSERT(fdata->getDepth() == cdata->getDepth());
#endif
    const int data_depth = fdata->getDepth();

    const Box<NDIM>& patch_box_fine = fine.getBox();
    const Index<NDIM>& patch_lower_fine = patch_box_fine.lower();
    Pointer<CartesianPatchGeometry<NDIM> > pgeom_fine = fine.getPatchGeometry();
    const double* const XLower_fine = pgeom_fine->getXLower();
    const double* const dx_fine = pgeom_fine->getDx();

    const Box<NDIM>& patch_box_crse = coarse.getBox();
    const Index<NDIM>& patch_lower_crse = patch_box_crse.lower();
    Pointer<CartesianPatchGeometry<NDIM> > pgeom_crse = coarse.getPatchGeometry();
    const double* const XLower_crse = pgeom_crse->getXLower();
    const double* const dx_crse = pgeom_crse->getDx();

    // Set all values in the fine box via quadratic interpolation from the
    // overlying coarse grid data.
    for (Box<NDIM>::Iterator b(fine_box); b; b++)
    {
        const Index<NDIM>& i_fine = b();
        const Index<NDIM> i_crse = coarsen(i_fine, ratio);

        // Determine the interpolation stencil in the coarse index space.
        Box<NDIM> stencil_box_crse(i_crse, i_crse);
        stencil_box_crse.grow(IntVector<NDIM>(1));

        // Determine the interpolation weights.
        static const int degree = 2;
        boost::array<boost::array<double, degree + 1>, NDIM> wgts(
            array_constant<boost::array<double, degree + 1>, NDIM>(
                boost::array<double, degree + 1>(array_constant<double, degree + 1>(0.0))));
        for (unsigned int axis = 0; axis < NDIM; ++axis)
        {
            const double X =
                XLower_fine[axis] + dx_fine[axis] * (static_cast<double>(i_fine(axis) - patch_lower_fine(axis)) + 0.5);
            std::vector<double> X_crse(degree + 1, 0.0);
            for (int i_crse = stencil_box_crse.lower()(axis), k = 0; i_crse <= stencil_box_crse.upper()(axis);
                 ++i_crse, ++k)
            {
                X_crse[k] =
                    XLower_crse[axis] + dx_crse[axis] * (static_cast<double>(i_crse - patch_lower_crse(axis)) + 0.5);
            }
            wgts[axis][0] = ((X - X_crse[1]) * (X - X_crse[2])) / ((X_crse[0] - X_crse[1]) * (X_crse[0] - X_crse[2]));
            wgts[axis][1] = ((X - X_crse[0]) * (X - X_crse[2])) / ((X_crse[1] - X_crse[0]) * (X_crse[1] - X_crse[2]));
            wgts[axis][2] = ((X - X_crse[0]) * (X - X_crse[1])) / ((X_crse[2] - X_crse[0]) * (X_crse[2] - X_crse[1]));
        }

        // Interpolate from the coarse grid to the fine grid.
        Index<NDIM> i_intrp;
        for (int d = 0; d < data_depth; ++d)
        {
            (*fdata)(i_fine, d) = 0.0;
#if (NDIM > 2)
            for (int i2 = 0; i2 <= degree; ++i2)
            {
                const double& wgt2 = wgts[2][i2];
                i_intrp(2) = stencil_box_crse.lower()(2) + i2;
#else
            const double wgt2 = 1.0;
#endif
#if (NDIM > 1)
                for (int i1 = 0; i1 <= degree; ++i1)
                {
                    const double& wgt1 = wgts[1][i1];
                    i_intrp(1) = stencil_box_crse.lower()(1) + i1;
#else
            const double wgt1 = 1.0;
#endif
                    for (int i0 = 0; i0 <= degree; ++i0)
                    {
                        const double& wgt0 = wgts[0][i0];
                        i_intrp(0) = stencil_box_crse.lower()(0) + i0;

                        (*fdata)(i_fine, d) += wgt0 * wgt1 * wgt2 * (*cdata)(i_intrp, d);
                    }
#if (NDIM > 1)
                }
#endif
#if (NDIM > 2)
            }
#endif
        }
    }
    return;
} // refine