fix_dof.h

#ifndef INCG_LF_FIXDOF_H
#define INCG_LF_FIXDOF_H
/***************************************************************************
 * LehrFEM++ - A simple C++ finite element libray for teaching
 * Developed from 2018 at the Seminar of Applied Mathematics of ETH Zurich,
 * lead developers Dr. R. Casagrande and Prof. R. Hiptmair
 ***************************************************************************/
 
#include "coomatrix.h"
 
namespace lf::assemble {
 
template <typename SCALAR, typename SELECTOR, typename RHSVECTOR>
void FixFlaggedSolutionComponents(SELECTOR &&selectvals, COOMatrix<SCALAR> &A,
                                  RHSVECTOR &b) {
  const lf::assemble::size_type N(A.cols());
  LF_ASSERT_MSG(A.rows() == N, "Matrix must be square!");
  LF_ASSERT_MSG(N == b.size(),
                "Mismatch N = " << N << " <-> b.size() = " << b.size());
  {
    // Multiply sparse matrix with the vector of fixed components and subtract
    // result from right hand side
    // To skip irrelevant components of fixed_vec rely on a temporary vector
    Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> tmp_vec(N);
    for (lf::assemble::gdof_idx_t k = 0; k < N; ++k) {
      const auto selval{selectvals(k)};
      if (selval.first) {
        tmp_vec[k] = selval.second;
      } else {
        tmp_vec[k] = SCALAR();
      }
    }
    A.MatVecMult(-1.0, tmp_vec, b);
  }
  // Set vector components of right-hand-side vector for prescribed values
  for (lf::assemble::gdof_idx_t k = 0; k < N; ++k) {
    const auto selval{selectvals(k)};
    if (selval.first) {
      b[k] = selval.second;
    }
  }
  // Set rows and columns of the sparse matrix corresponding to the fixed
  // solution components to zero
  A.setZero([&selectvals](gdof_idx_t i, gdof_idx_t j) {
    return (selectvals(i).first || selectvals(j).first);
  });
  // Old implementation, demonstrating what is going on
  // lf::assemble::COOMatrix<double>::TripletVec::iterator new_last =
  //     std::remove_if(
  //         A.triplets().begin(), A.triplets().end(),
  //         [&fixed_comp_flags](
  //             typename lf::assemble::COOMatrix<double>::Triplet &triplet) {
  //           return (fixed_comp_flags[triplet.row()] ||
  //                   fixed_comp_flags[triplet.col()]);
  //         });
  // // Adjust size of triplet vector
  // A.triplets().erase(new_last, A.triplets().end());
  // Add Unit diagonal entries corrresponding to fixed components
  for (lf::assemble::gdof_idx_t dofnum = 0; dofnum < N; ++dofnum) {
    const auto selval{selectvals(dofnum)};
    if (selval.first) {
      A.AddToEntry(dofnum, dofnum, 1.0);
    }
  }
}
 
template <typename SCALAR, typename SELECTOR, typename RHSVECTOR>
void FixFlaggedSolutionCompAlt(SELECTOR &&selectvals, COOMatrix<SCALAR> &A,
                               RHSVECTOR &b) {
  const lf::assemble::size_type N(A.cols());
  LF_ASSERT_MSG(A.rows() == N, "Matrix must be square!");
  LF_ASSERT_MSG(N == b.size(),
                "Mismatch N = " << N << " <-> b.size() = " << b.size());
 
  // Set vector components of right-hand-side vector for prescribed values
  for (lf::assemble::gdof_idx_t k = 0; k < N; ++k) {
    const auto selval{selectvals(k)};
    if (selval.first) {
      b[k] = selval.second;
    }
  }
  // Set rows and columns of the sparse matrix corresponding to the fixed
  // solution components to zero
  A.setZero([&selectvals](gdof_idx_t i, gdof_idx_t /*unused*/) {
    return (selectvals(i).first);
  });
  // Old implementation showing the algorithm:
  // lf::assemble::COOMatrix<double>::TripletVec::iterator new_last =
  //     std::remove_if(
  //         A.triplets().begin(), A.triplets().end(),
  //         [&fixed_comp_flags](
  //             typename lf::assemble::COOMatrix<double>::Triplet &triplet) {
  //           return (fixed_comp_flags[triplet.row()]);
  //         });
  // // Adjust size of triplet vector
  // A.triplets().erase(new_last, A.triplets().end());
  // Add Unit diagonal entries corrresponding to fixed components
  for (lf::assemble::gdof_idx_t dofnum = 0; dofnum < N; ++dofnum) {
    if (selectvals(dofnum).first) {
      A.AddToEntry(dofnum, dofnum, 1.0);
    }
  }
}
 
template <typename SCALAR>
using fixed_components_t =
    std::vector<std::pair<lf::assemble::gdof_idx_t, SCALAR>>;
 
template <typename SCALAR, typename RHSVECTOR>
void FixSolutionComponentsLse(
    const fixed_components_t<SCALAR> &fixed_components, COOMatrix<SCALAR> &A,
    RHSVECTOR &b) {
  const lf::assemble::size_type N(A.cols());
  LF_ASSERT_MSG(A.rows() == N, "Matrix must be square!");
  LF_ASSERT_MSG(N == b.size(),
                "Mismatch N = " << N << " <-> b.size() = " << b.size());
  // Set up a vector containing the negative fixed solution components and a
  // flag vector, whose entry is `true` if the corresponding
  // vector component is meant to be fixed
  Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> fixed_vec(b.size());
  fixed_vec.setZero();
  std::vector<bool> fixed_comp_flags(N, false);
  for (typename fixed_components_t<SCALAR>::const_reference idx_val_pair :
       fixed_components) {
    LF_ASSERT_MSG(idx_val_pair.first < N,
                  "Index " << idx_val_pair.first << " >= N = " << N);
    fixed_vec[idx_val_pair.first] += idx_val_pair.second;
    fixed_comp_flags[idx_val_pair.first] = true;
  }
  // FixFlaggedSolutionComponents<SCALAR>(fixed_comp_flags, fixed_vec, A, b);
  FixFlaggedSolutionCompAlt<SCALAR>(
      [&fixed_comp_flags,
       &fixed_vec](lf::assemble::gdof_idx_t i) -> std::pair<bool, double> {
        LF_ASSERT_MSG((i < fixed_comp_flags.size()) && (i < fixed_vec.size()),
                      "Illegal index " << i);
        return std::make_pair(fixed_comp_flags[i], fixed_vec[i]);
      },
      A, b);
}
 
}  // namespace lf::assemble
 
#endif