Method calculateVolumeElement

src/IB/ConstraintIBMethod.cpp:1383–1528

} // calculateMomentumOfKinematicsVelocity

void
ConstraintIBMethod::calculateVolumeElement()
{
    using StructureParameters = ConstraintIBKinematics::StructureParameters;

    // Initialize variables and variable contexts associated with Eulerian
    // tracking of the Lagrangian points.
    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
    const IntVector<NDIM> cell_ghosts = 0;
    Pointer<CellVariable<NDIM, int>> vol_cc_var = new CellVariable<NDIM, int>(d_object_name + "::vol_cc_var");
    const int vol_cc_scratch_idx = var_db->registerVariableAndContext(vol_cc_var, d_scratch_context, cell_ghosts);

    const int coarsest_ln = 0;
    const int finest_ln = d_hierarchy->getFinestLevelNumber();
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;
        Pointer<PatchLevel<NDIM>> level = d_hierarchy->getPatchLevel(ln);
        level->allocatePatchData(vol_cc_scratch_idx, 0.0);
        for (PatchLevel<NDIM>::Iterator p(level); p; p++)
        {
            Pointer<Patch<NDIM>> patch = level->getPatch(p());
            Pointer<CellData<NDIM, int>> vol_cc_scratch_idx_data = patch->getPatchData(vol_cc_scratch_idx);
            vol_cc_scratch_idx_data->fill(0, 0);
        }
    }

    const int lag_node_index_idx = d_l_data_manager->getLNodePatchDescriptorIndex();
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;

        // Get LData corresponding to the present position of the structures.
        const boost::multi_array_ref<double, 2>& X_data = *d_l_data_manager->getLData("X", ln)->getLocalFormVecArray();
        Pointer<PatchLevel<NDIM>> level = d_hierarchy->getPatchLevel(ln);

        // Get structures on this level.
        const std::vector<int> structIDs = d_l_data_manager->getLagrangianStructureIDs(ln);
        const size_t structs_on_this_ln = structIDs.size();
        for (size_t struct_no = 0; struct_no < structs_on_this_ln; ++struct_no)
        {
            // If the volume element has already been set, then skip the volume
            // computation
            if (d_vol_element_is_set[struct_no])
            {
                tbox::plog << "Skipping volume element computation for structure no. " << struct_no << std::endl;
                tbox::pout << "Skipping volume element computation for structure no. " << struct_no << std::endl;
                continue;
            }
            std::pair<int, int> lag_idx_range =
                d_l_data_manager->getLagrangianStructureIndexRange(structIDs[struct_no], ln);
            Pointer<ConstraintIBKinematics> ptr_ib_kinematics =
                *std::find_if(d_ib_kinematics.begin(), d_ib_kinematics.end(), find_struct_handle(lag_idx_range));
            const int location_struct_handle =
                find_struct_handle_position(d_ib_kinematics.begin(), d_ib_kinematics.end(), ptr_ib_kinematics);

            for (PatchLevel<NDIM>::Iterator p(level); p; p++)
            {