Method computeL2Projection

ibtk/src/lagrangian/FEProjector.cpp:676–770

}

bool
FEProjector::computeL2Projection(PetscVector<double>& U_vec,
                                 PetscVector<double>& F_vec,
                                 const std::string& system_name,
                                 const bool consistent_mass_matrix,
                                 const bool close_U,
                                 const bool close_F,
                                 const double tol,
                                 const unsigned int max_its)
{
    tbox::Pointer<tbox::Timer>& system_timer = d_linear_solve_system_timers[system_name];
    if (system_timer.isNull())
    {
        system_timer =
            TimerManager::getManager()->getTimer("IBTK::FEProjector::computeL2Projection()[" + system_name + "]");
        TBOX_ASSERT(system_timer);
    }
    IBTK_TIMER_START(system_timer);

    int ierr;
    bool converged = false;

    if (close_F) F_vec.close();
    const System& system = d_fe_data->getEquationSystems()->get_system(system_name);

    const MeshBase& mesh = d_fe_data->getEquationSystems()->get_mesh();
    const unsigned int dim = mesh.mesh_dimension();

    FEType fe_type = system.get_dof_map().variable_type(0);

    // We can use the diagonal mass matrix directly if we do not need a
    // consistent mass matrix *and* there are no constraints.
    //
    // TODO: this would also work with Dirichlet boundary constraints but it is not
    // as easy to detect those constraints.
    if (!consistent_mass_matrix && system.get_dof_map().n_constrained_dofs() == 0)
    {
        PetscVector<double>* M_diag_vec = buildDiagonalL2MassMatrix(system_name);
        ierr = VecPointwiseDivide(U_vec.vec(), F_vec.vec(), M_diag_vec->vec());
        IBTK_CHKERRQ(ierr);
        converged = true;
    }
    else
    {
        std::pair<PetscLinearSolver<double>*, PetscMatrix<double>*> proj_solver_components =
            consistent_mass_matrix ? buildL2ProjectionSolver(system_name) : buildLumpedL2ProjectionSolver(system_name);
        PetscMatrix<double>& lumped_mass = *buildLumpedL2ProjectionSolver(system_name).second;
        PetscLinearSolver<double>* solver = proj_solver_components.first;
        PetscMatrix<double>* M_mat = proj_solver_components.second;
        PetscBool rtol_set;
        double runtime_rtol;
        ierr = PetscOptionsGetReal(nullptr, "", "-ksp_rtol", &runtime_rtol, &rtol_set);
        IBTK_CHKERRQ(ierr);
        PetscBool max_it_set;
        int runtime_max_it;
        ierr = PetscOptionsGetInt(nullptr, "", "-ksp_max_it", &runtime_max_it, &max_it_set);
        IBTK_CHKERRQ(ierr);
        ierr = KSPSetFromOptions(solver->ksp());