Method computeForcing

src/IB/BrinkmanAdvDiffBcHelper.cpp:629–881

} // computeDiffusionCoefficient

void
BrinkmanAdvDiffBcHelper::computeForcing(int F_idx, Pointer<CellVariable<NDIM, double>> Q_var)
{
#if !defined(NDEBUG)
    TBOX_ASSERT(d_Q_bc.find(Q_var) != d_Q_bc.end());
#endif
    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
    std::vector<BCProperties> brinkman_zones = d_Q_bc[Q_var];

    Pointer<PatchHierarchy<NDIM>> patch_hierarchy = d_adv_diff_solver->getPatchHierarchy();
    const int coarsest_ln = 0;
    const int finest_ln = patch_hierarchy->getFinestLevelNumber();

    HierarchyCellDataOpsReal<NDIM, double> hier_cc_data_ops(patch_hierarchy, coarsest_ln, finest_ln);
    hier_cc_data_ops.setToScalar(F_idx, 0.0);

    // First, deal with Dirichlet, spatially constant g, which can be computed pointwise.
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        Pointer<PatchLevel<NDIM>> level = patch_hierarchy->getPatchLevel(ln);
        for (PatchLevel<NDIM>::Iterator p(level); p; p++)
        {
            Pointer<Patch<NDIM>> patch = level->getPatch(p());
            const Pointer<CartesianPatchGeometry<NDIM>> patch_geom = patch->getPatchGeometry();
            const Box<NDIM>& patch_box = patch->getBox();
            const double* patch_dx = patch_geom->getDx();
            double vol_cell = 1.0;
            for (int d = 0; d < NDIM; ++d) vol_cell *= patch_dx[d];

            Pointer<CellData<NDIM, double>> F_data = patch->getPatchData(F_idx);

            for (Box<NDIM>::Iterator it(patch_box); it; it++)
            {
                CellIndex<NDIM> ci(it());

                // Loop over all the level sets and add up their contributions.
                // TODO: Pointwise check that level sets don't overlap
                double brinkman_forcing = 0.0;
                for (const auto& bc_prop : brinkman_zones)
                {
                    // Get the BC specifications for each zone.
                    Pointer<CellVariable<NDIM, double>> ls_solid_var = bc_prop.ls_solid_var;
                    AdvDiffBrinkmanPenalizationBcType bc_type = bc_prop.bc_type;
                    double bc_val = bc_prop.bc_val;
                    double eta = bc_prop.eta;

                    bool requires_callback = bc_prop.callback != nullptr;
                    double num_interface_cells = bc_prop.num_interface_cells;
                    const double alpha = num_interface_cells * std::pow(vol_cell, 1.0 / static_cast<double>(NDIM));

                    if (requires_callback) continue;
#if !defined(NDEBUG)
                    TBOX_ASSERT(!requires_callback);
#endif
                    const int phi_idx =
                        var_db->mapVariableAndContextToIndex(ls_solid_var, d_adv_diff_solver->getNewContext());
                    Pointer<CellData<NDIM, double>> ls_solid_data = patch->getPatchData(phi_idx);
                    double phi = (*ls_solid_data)(ci);