Method initializeMassDataOnPatchLevel

src/IB/IBRedundantInitializer.cpp:1112–1170

} // initializeDataOnPatchLevel

unsigned int
IBRedundantInitializer::initializeMassDataOnPatchLevel(const unsigned int /*global_index_offset*/,
                                                       const unsigned int local_index_offset,
                                                       Pointer<LData> M_data,
                                                       Pointer<LData> K_data,
                                                       const Pointer<PatchHierarchy<NDIM>> hierarchy,
                                                       const int level_number,
                                                       const double /*init_data_time*/,
                                                       const bool /*can_be_refined*/,
                                                       const bool /*initial_time*/,
                                                       LDataManager* const /*l_data_manager*/)
{
#if !defined(NDEBUG)
    TBOX_ASSERT(d_data_processed);
#endif

    // Loop over all patches in the specified level of the patch level and
    // initialize the local vertices.
    boost::multi_array_ref<double, 1>& M_array = *M_data->getLocalFormArray();
    boost::multi_array_ref<double, 1>& K_array = *K_data->getLocalFormArray();
    int local_idx = invalid_index;
    int local_node_count = 0;
    Pointer<PatchLevel<NDIM>> level = hierarchy->getPatchLevel(level_number);
    for (PatchLevel<NDIM>::Iterator p(level); p; p++)
    {
        Pointer<Patch<NDIM>> patch = level->getPatch(p());

        // Initialize the vertices whose initial locations will be within the
        // given patch.
        std::vector<std::pair<int, int>> patch_vertices;
        getPatchVertices(patch_vertices, patch, hierarchy);
        local_node_count += patch_vertices.size();
        for (const auto& point_idx : patch_vertices)
        {
            const int local_petsc_idx = ++local_idx + local_index_offset;

            // Initialize the mass and penalty stiffness coefficient
            // corresponding to the present vertex.
            const BdryMassSpec& spec = getVertexBdryMassSpec(point_idx, level_number);
            const double M = spec.bdry_mass;
            const double K = spec.stiffness;

            // Avoid division by zero at massless nodes.
            if (IBTK::abs_equal_eps(M, 0.0))
            {
                M_array[local_petsc_idx] = std::numeric_limits<double>::epsilon();
                K_array[local_petsc_idx] = 0.0;
            }
            else
            {
                M_array[local_petsc_idx] = M;
                K_array[local_petsc_idx] = K;
            }
        }
    }
    M_data->restoreArrays();
    K_data->restoreArrays();
    return local_node_count;