Method spreadFluidSource

src/IB/IBMethod.cpp:1068–1223

} // computeLagrangianFluidSource

void
IBMethod::spreadFluidSource(const int q_data_idx,
                            RobinPhysBdryPatchStrategy* /*q_phys_bdry_op*/,
                            const std::vector<Pointer<RefineSchedule<NDIM>>>& /*q_prolongation_scheds*/,
                            const double data_time)
{
    if (!d_ib_source_fcn) return;

    const int coarsest_ln = 0;
    const int finest_ln = d_hierarchy->getFinestLevelNumber();

    // Get the present source locations.
    std::vector<Pointer<LData>>* X_data;
    bool* X_needs_ghost_fill;
    getLECouplingPositionData(&X_data, &X_needs_ghost_fill, data_time);
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        if (d_n_src[ln] == 0) continue;
        d_ib_source_fcn->getSourceLocations(
            d_X_src[ln], d_r_src[ln], (*X_data)[ln], d_hierarchy, ln, data_time, d_l_data_manager);
    }

    // Spread the sources/sinks onto the Cartesian grid.
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        if (d_n_src[ln] == 0) continue;
#if !defined(NDEBUG)
        TBOX_ASSERT(ln == d_hierarchy->getFinestLevelNumber());
#endif
        Pointer<PatchLevel<NDIM>> level = d_hierarchy->getPatchLevel(ln);
        const IntVector<NDIM>& ratio = level->getRatio();
        const Pointer<CartesianGridGeometry<NDIM>> grid_geom = level->getGridGeometry();
        for (PatchLevel<NDIM>::Iterator p(level); p; p++)
        {
            Pointer<Patch<NDIM>> patch = level->getPatch(p());
            const Box<NDIM>& patch_box = patch->getBox();
            const hier::Index<NDIM>& patch_lower = patch_box.lower();
            const Pointer<CartesianPatchGeometry<NDIM>> pgeom = patch->getPatchGeometry();
            const double* const xLower = pgeom->getXLower();
            const double* const dx = pgeom->getDx();
            const Pointer<CellData<NDIM, double>> q_data = patch->getPatchData(q_data_idx);
            for (int n = 0; n < d_n_src[ln]; ++n)
            {
                // The source radius must be an integer multiple of the grid
                // spacing.
                std::array<double, NDIM> r;
                for (unsigned int d = 0; d < NDIM; ++d)
                {
                    r[d] = std::max(std::floor(d_r_src[ln][n] / dx[d] + 0.5), 2.0) * dx[d];
                }

                // Determine the approximate source stencil box.
                const hier::Index<NDIM> i_center = IndexUtilities::getCellIndex(d_X_src[ln][n], grid_geom, ratio);
                Box<NDIM> stencil_box(i_center, i_center);
                for (unsigned int d = 0; d < NDIM; ++d)
                {
                    stencil_box.grow(d, static_cast<int>(r[d] / dx[d]) + 1);
                }