Method solveSystem

src/IB/CIBStaggeredStokesSolver.cpp:256–377

} // setPhysicalBoundaryHelper

bool
CIBStaggeredStokesSolver::solveSystem(SAMRAIVectorReal<NDIM, double>& x, SAMRAIVectorReal<NDIM, double>& b)
{
    // Create packaged vectors for the Saddle point solver.
    d_x_wide->copyVector(Pointer<SAMRAIVectorReal<NDIM, double>>(&x, false));
    d_b_wide->copyVector(Pointer<SAMRAIVectorReal<NDIM, double>>(&b, false));

    // Wrap SAMRAI vector into PETSc Vec
    Vec u_p = PETScSAMRAIVectorReal::createPETScVector(d_x_wide);
    Vec g_h = PETScSAMRAIVectorReal::createPETScVector(d_b_wide);

    // Get the Lagrange multiplier that maintains the rigidity constraint.
    // NOTE: We need L at new time to solve for it.
    Vec L;
    d_cib_strategy->getConstraintForce(&L, d_new_time);

    // Get the free body velocities to solve for.
    // NOTE: We need U at new time in the solver.
    Vec U;
    d_cib_strategy->getFreeRigidVelocities(&U, d_new_time);

    // Set the imposed velocity for all bodies in the RHS.
    Vec V;
    VecDuplicate(L, &V);
    for (unsigned part = 0; part < d_num_rigid_parts; ++part)
    {
        RigidDOFVector U_part;
        d_cib_strategy->getNewRigidBodyVelocity(part, U_part);

        // Zero-out free velocities.
        int num_free_dofs = 0;
        const FreeRigidDOFVector& solve_dofs = d_cib_strategy->getSolveRigidBodyVelocity(part, num_free_dofs);
        for (int k = 0; k < s_max_free_dofs; ++k)
        {
            if (solve_dofs[k]) U_part[k] = 0.0;
        }

        const double interp_scale = d_sp_solver->getInterpScale();
        U_part *= -interp_scale;
        d_cib_strategy->setRigidBodyVelocity(part, U_part, V);
    }

    // Get the net external force and torque on the bodies.
    Vec F;
    d_cib_strategy->getNetExternalForceTorque(&F, d_new_time);

    // Create multivector to pass it to the saddle point solver.
    std::vector<Vec> vx(3), vb(3);
    vx[0] = u_p;
    vx[1] = L;
    vx[2] = U;
    vb[0] = g_h;
    vb[1] = V;
    vb[2] = F;

    Vec mv_x, mv_b;
    VecCreateNest(PETSC_COMM_WORLD, 3, nullptr, &vx[0], &mv_x);
    VecCreateNest(PETSC_COMM_WORLD, 3, nullptr, &vb[0], &mv_b);