Method updateStructurePositionMidPointStep

src/IB/ConstraintIBMethod.cpp:2228–2321

} // eulerStep

void
ConstraintIBMethod::updateStructurePositionMidPointStep()
{
    using StructureParameters = ConstraintIBKinematics::StructureParameters;
    const int coarsest_ln = 0;
    const int finest_ln = d_hierarchy->getFinestLevelNumber();
    const double dt = d_FuRMoRP_new_time - d_FuRMoRP_current_time;

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

        boost::multi_array_ref<double, 2>& X_new_MidPoint_data = *d_l_data_X_new_MidPoint[ln]->getLocalFormVecArray();
        const boost::multi_array_ref<double, 2>& X_current_data =
            *d_l_data_manager->getLData("X", ln)->getLocalFormVecArray();
        const boost::multi_array_ref<double, 2>& U_half_data = *d_l_data_U_half[ln]->getLocalFormVecArray();

        const Pointer<LMesh> mesh = d_l_data_manager->getLMesh(ln);
        const std::vector<LNode*>& local_nodes = mesh->getLocalNodes();

        // 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)
        {
            std::pair<int, int> lag_idx_range =
                d_l_data_manager->getLagrangianStructureIndexRange(structIDs[struct_no], ln);
            const int offset = lag_idx_range.first;
            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);
            const StructureParameters& struct_param = ptr_ib_kinematics->getStructureParameters();
            const std::string position_update_method = struct_param.getPositionUpdateMethod();
            const std::vector<std::vector<double>>& new_shape = ptr_ib_kinematics->getShape(ln);

            for (const auto& node_idx : local_nodes)
            {
                const int lag_idx = node_idx->getLagrangianIndex();
                if (lag_idx_range.first <= lag_idx && lag_idx < lag_idx_range.second)
                {
                    const int local_idx = node_idx->getLocalPETScIndex();
                    const double* const U_half = &U_half_data[local_idx][0];
                    const double* const X_current = &X_current_data[local_idx][0];
                    double* const X_new = &X_new_MidPoint_data[local_idx][0];

                    if (position_update_method == "CONSTRAINT_VELOCITY")
                    {
                        for (int d = 0; d < NDIM; ++d)
                        {
                            X_new[d] = X_current[d] + dt * U_half[d];
                        }
                    }
                    else if (position_update_method == "CONSTRAINT_POSITION")
                    {
                        for (int d = 0; d < NDIM; ++d)