Assemble Schur complement (Implementation) More...

Inheritance diagram for MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >:

Collaboration diagram for MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >:

Public Types
using	OP = OP_SCHUR_ASSEMBLE_BASE

Public Member Functions
	OpSchurAssembleEndImpl (std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > schur_ao, SmartPetscObj< Mat > schur_mat, bool sym_schur, bool symm_op)
	Construct a new Op Schur Assemble End object.

Protected Member Functions
template<typename I >
MoFEMErrorCode	doWorkImpl (int side, EntityType type, EntitiesFieldData::EntData &data)

Protected Attributes
std::vector< std::string >	fieldsName

std::vector< boost::shared_ptr< Range > >	fieldEnts

SmartPetscObj< AO >	schurAO

SmartPetscObj< Mat >	schurMat

bool	symSchur

MatrixDouble	blockMat

MatrixDouble	invMat

MatrixDouble	bM

MatrixDouble	abM

MatrixDouble	abcM

Detailed Description

template<typename OP_SCHUR_ASSEMBLE_BASE>
struct MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >

Assemble Schur complement (Implementation)

Definition at line 88 of file Schur.cpp.

Member Typedef Documentation

◆ OP

template<typename OP_SCHUR_ASSEMBLE_BASE >

using MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >::OP = OP_SCHUR_ASSEMBLE_BASE

Definition at line 90 of file Schur.cpp.

Constructor & Destructor Documentation

◆ OpSchurAssembleEndImpl()

template<typename OP_SCHUR_ASSEMBLE_BASE >

MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >::OpSchurAssembleEndImpl	(	std::vector< std::string >	fields_name,
		std::vector< boost::shared_ptr< Range > >	field_ents,
		SmartPetscObj< AO >	schur_ao,
		SmartPetscObj< Mat >	schur_mat,
		bool	sym_schur,
		bool	symm_op
	)

Construct a new Op Schur Assemble End object.

Parameters

fields_name	list of fields
field_ents	list of entities on which schur complement is applied (can be empty)
schur_ao	map indices to schur matrix indices
schur_mat	schur matrix
sym_schur	true if schur is symmetric
symm_op	true if block diagonal is symmetric

Definition at line 583 of file Schur.cpp.

    : OP(NOSPACE, OP::OPSPACE, symm_op), fieldsName(fields_name),
      fieldEnts(field_ents), schurAO(schur_ao), schurMat(schur_mat),
      symSchur(sym_schur) {}

Member Function Documentation

◆ doWorkImpl()

template<typename OP_SCHUR_ASSEMBLE_BASE >

template<typename I >

MoFEMErrorCode MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >::doWorkImpl	(	int	side,
		EntityType	type,
		EntitiesFieldData::EntData &	data
	)

protected

Definition at line 594 of file Schur.cpp.

                                                               {
 
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  PetscLogEventBegin(SchurEvents::MOFEM_EVENT_opSchurAssembleEnd, 0, 0, 0, 0);
#endif
 
#ifndef NDEBUG
  if constexpr (debug_schur)
    MOFEM_LOG("SELF", Sev::noisy) << "Schur assemble begin -> end";
#endif
 
  auto get_a00_uids = [&]() {
    auto get_field_bit = [&](auto &name) {
      return OP::getPtrFE()->mField.get_field_bit_number(name);
    };
 
    std::vector<std::pair<UId, UId>> a00_uids;
    a00_uids.reserve(fieldsName.size());
    for (auto ss = 0; ss != fieldsName.size(); ++ss) {
      auto field_bit = get_field_bit(fieldsName[ss]);
      auto row_ents = fieldEnts[ss];
      if (row_ents) {
        for (auto p = row_ents->pair_begin(); p != row_ents->pair_end(); ++p) {
          auto lo_uid =
              FieldEntity::getLoLocalEntityBitNumber(field_bit, p->first);
          auto hi_uid =
              FieldEntity::getHiLocalEntityBitNumber(field_bit, p->second);
          a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
        }
      } else {
        auto lo_uid = FieldEntity::getLoLocalEntityBitNumber(
            field_bit, get_id_for_min_type<MBVERTEX>());
        auto hi_uid = FieldEntity::getHiLocalEntityBitNumber(
            field_bit, get_id_for_max_type<MBENTITYSET>());
        a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
      }
    }
    return a00_uids;
  };
 
#ifndef NDEBUG
  auto get_field_name = [&](auto uid) {
    return OP::getPtrFE()->mField.get_field_name(field_bit_from_bit_number(
        FieldEntity::getFieldBitNumberFromUniqueId(uid)));
  };
  auto list_storage = [&](auto &storage) {
    MoFEMFunctionBegin;
    int i = 0;
    for (auto &p : storage) {
      MOFEM_LOG("SELF", Sev::noisy)
          << "List schur storage: " << i << " " << p->iDX << ": "
          << get_field_name(p->uidRow) << " " << get_field_name(p->uidCol)
          << " : " << p->getMat().size1() << " " << p->getMat().size2();
      ++i;
    }
    MoFEMFunctionReturn(0);
  };
#endif // NDEBUG
 
  auto assemble_dense_blocks = [&]() {
    using matrix_range = ublas::matrix_range<MatrixDouble>;
    using range = ublas::range;
    MoFEMFunctionBegin;
    auto &storage = SchurElemMats::schurL2Storage;
 
    auto assemble_schur = [this](auto &m, auto &uid_row, auto &uid_col,
                                 auto *row_ind_ptr, auto *col_ind_ptr) {
      MoFEMFunctionBegin;
 
      if (schurMat) {
 
        if (auto ierr = this->assembleSchurMat(
 
                schurMat, uid_row, *row_ind_ptr, uid_col, *col_ind_ptr, m,
                ADD_VALUES
 
                )) {
#ifndef NDEBUG
          auto field_ents = OP::getPtrFE()->mField.get_field_ents();
          auto row_ent_it = field_ents->find(uid_row);
          auto col_ent_it = field_ents->find(uid_col);
          MOFEM_LOG_CHANNEL("SELF");
          if (row_ent_it != field_ents->end())
            MOFEM_LOG("SELF", Sev::error)
                << "Assemble row entity: " << (*row_ent_it)->getName() << " "
                << (*col_ent_it)->getEntTypeName() << " side "
                << (*row_ent_it)->getSideNumber();
          if (col_ent_it != field_ents->end())
            MOFEM_LOG("SELF", Sev::error)
                << "Assemble col entity: " << (*col_ent_it)->getName() << " "
                << (*col_ent_it)->getEntTypeName() << " side "
                << (*col_ent_it)->getSideNumber();
#endif // NDEBUG
          CHK_THROW_MESSAGE(ierr, "MatSetValues");
        }
      }
 
      MoFEMFunctionReturn(0);
    };
 
    auto a00_uids = get_a00_uids();
 
    auto get_block_indexing = [&](auto &a00_uids) {
      // iterate over a00 uids and find blocks
      std::vector<const SchurElemMats *> block_list;
      block_list.reserve(storage.size());
      for (auto &rp_uid : a00_uids) {
        auto [rlo_uid, rhi_uid] = rp_uid;
        for (auto &cp_uid : a00_uids) {
          auto [clo_uid, chi_uid] = cp_uid;
 
          auto it =
              storage.template get<SchurElemMats::uid_mi_tag>().lower_bound(
                  boost::make_tuple(rlo_uid, clo_uid));
          auto hi_it =
              storage.template get<SchurElemMats::uid_mi_tag>().upper_bound(
                  boost::make_tuple(rhi_uid, chi_uid));
 
          for (; it != hi_it; ++it) {
            if ((*it)->uidRow >= rlo_uid && (*it)->uidRow < rhi_uid &&
                (*it)->uidCol >= clo_uid && (*it)->uidCol < chi_uid) {
              block_list.push_back(*it);
            }
          }
        }
      }
 
      // create block indexes map for blockMat
      std::map<UId, std::pair<size_t, const VectorInt *>>
          block_indexing; // uid block map
      for (auto d : block_list) {
        if (block_indexing.find(d->uidRow) == block_indexing.end()) {
          block_indexing[d->uidRow] =
              std::make_pair(d->getRowInd().size(), &(d->getRowInd()));
        }
      }
 
      // set indexes to block
      int mat_block_size = 0; // size of block matrix
      for (auto &p_uid : a00_uids) {
        auto [lo_uid, hi_uid] = p_uid;
        auto lo = block_indexing.lower_bound(lo_uid);
        auto up = block_indexing.upper_bound(hi_uid);
        for (; lo != up; ++lo) {
          lo->second.first = mat_block_size;
          mat_block_size += lo->second.second->size();
        }
      }
 
      return std::make_tuple(block_list, block_indexing, mat_block_size);
    };
 
    auto get_schur_block_list = [&](auto &block_indexing) {
      std::vector<const SchurElemMats *> block_list;
      block_list.reserve(storage.size());
      for (auto &s : storage) {
        if (block_indexing.find(s->uidRow) == block_indexing.end() &&
            block_indexing.find(s->uidCol) == block_indexing.end()) {
          block_list.push_back(s);
        }
      }
      return block_list;
    };
 
    auto bi = get_block_indexing(a00_uids);
    auto &[block_list, block_indexing, block_mat_size] = bi;
 
    if (block_mat_size == 0) {
 
      if (schurAO) {
        for (auto &s : storage) {
          CHKERR AOApplicationToPetsc(schurAO, s->getRowInd().size(),
                                      &*s->getRowInd().begin());
          CHKERR AOApplicationToPetsc(schurAO, s->getColInd().size(),
                                      &*s->getColInd().begin());
        }
      }
 
      for (auto &s : storage) {
        auto &m = s->getMat();
        CHKERR assemble_schur(m, s->uidRow, s->uidCol, &(s->getRowInd()),
                              &(s->getColInd()));
      }
      MoFEMFunctionReturnHot(0);
    }
 
    blockMat.resize(block_mat_size, block_mat_size, false);
    blockMat.clear();
 
    auto get_range = [](auto &bi) {
      return range(bi.first, bi.first + bi.second->size());
    };
 
    for (auto &s : block_list) {
      auto &m = s->getMat();
#ifndef NDEBUG
      if (block_indexing.find(s->uidRow) == block_indexing.end())
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong rlo_uid");
      if (block_indexing.find(s->uidCol) == block_indexing.end())
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong clo_uid");
#endif // NDEBUG
 
      auto &rbi = block_indexing.at(s->uidRow);
      auto &cbi = block_indexing.at(s->uidCol);
 
      auto sub_mat = matrix_range(blockMat, get_range(rbi), get_range(cbi));
      sub_mat = m;
    }
 
    auto get_zeroed_indices = [&](auto extractor_uid, auto extractor_ind) {
      auto &[block_list, block_indexing, block_mat_size] = bi;
      std::vector<int> zeroed_indices;
      zeroed_indices.reserve(block_mat_size);
      for (auto &s : block_list) {
        auto &bi = block_indexing.at(extractor_uid(s));
        auto &ind = extractor_ind(s);
        for (auto it = ind.begin(); it != ind.end(); ++it) {
          if (*it < 0) {
            auto idx = bi.first + std::distance(ind.begin(), it);
            zeroed_indices.push_back(idx);
          }
        }
      }
      std::sort(zeroed_indices.begin(), zeroed_indices.end());
      auto it = std::unique(zeroed_indices.begin(), zeroed_indices.end());
      zeroed_indices.resize(std::distance(zeroed_indices.begin(), it));
      return zeroed_indices;
    };
    auto zero_rows = get_zeroed_indices(
        [](auto &s) { return s->uidRow; },
        [](auto &s) -> VectorInt & { return s->getRowInd(); });
    auto zero_cols = get_zeroed_indices(
        [](auto &s) { return s->uidCol; },
        [](auto &s) -> VectorInt & { return s->getColInd(); });
 
    for (auto i : zero_rows) {
      for (auto j = 0; j != blockMat.size2(); ++j) {
        blockMat(i, j) = 0;
      }
    }
    for (auto j : zero_cols) {
      for (auto i = 0; i != blockMat.size1(); ++i) {
        blockMat(i, j) = 0;
      }
    }
    for (auto i : zero_rows) {
      blockMat(i, i) = 1;
    }
 
    CHKERR I::invertMat(blockMat, invMat);
 
    // clear storage and block list from a00 blocks, no more needed
    for (auto &s : block_list) {
      auto it = storage.template get<SchurElemMats::uid_mi_tag>().find(
          boost::make_tuple(s->uidRow, s->uidCol));
      storage.template get<SchurElemMats::uid_mi_tag>().erase(it);
    }
    block_list.clear();
 
    if (schurAO) {
      for (auto &s : storage) {
        if (block_indexing.find(s->uidRow) == block_indexing.end()) {
          CHKERR AOApplicationToPetsc(schurAO, s->getRowInd().size(),
                                      &*s->getRowInd().begin());
        }
        if (block_indexing.find(s->uidCol) == block_indexing.end()) {
          CHKERR AOApplicationToPetsc(schurAO, s->getColInd().size(),
                                      &*s->getColInd().begin());
        }
      }
    }
 
    auto schur_block_list = get_schur_block_list(block_indexing);
    for (auto &s : schur_block_list) {
      auto &m = s->getMat();
      CHKERR assemble_schur(m, s->uidRow, s->uidCol, &(s->getRowInd()),
                            &(s->getColInd()));
    }
    for (auto &s : schur_block_list) {
      auto it = storage.template get<SchurElemMats::uid_mi_tag>().find(
          boost::make_tuple(s->uidRow, s->uidCol));
      storage.template get<SchurElemMats::uid_mi_tag>().erase(it);
    }
    schur_block_list.clear();
 
    for (
 
        auto rp_uid_it = a00_uids.begin(); rp_uid_it != a00_uids.end();
        ++rp_uid_it
 
    ) {
      auto [rlo_uid, rhi_uid] = *rp_uid_it;
      for (auto rm = block_indexing.lower_bound(rlo_uid);
           rm != block_indexing.upper_bound(rhi_uid); ++rm) {
        auto &rbi = rm->second;
 
        auto a_lo_tmp =
            storage.template get<SchurElemMats::col_mi_tag>().lower_bound(
                rm->first);
        auto a_hi =
            storage.template get<SchurElemMats::col_mi_tag>().upper_bound(
                rm->first);
 
        for (
 
            auto cp_uid_it = (symSchur) ? rp_uid_it : a00_uids.begin();
            cp_uid_it != a00_uids.end(); ++cp_uid_it
 
        ) {
          auto [clo_uid, chi_uid] = *cp_uid_it;
          for (auto cm = block_indexing.lower_bound(clo_uid);
               cm != block_indexing.upper_bound(chi_uid); ++cm) {
            auto &cbi = cm->second;
 
            auto c_lo_tmp =
                storage.template get<SchurElemMats::uid_mi_tag>().lower_bound(
                    boost::make_tuple(cm->first, 0));
            auto c_hi =
                storage.template get<SchurElemMats::uid_mi_tag>().upper_bound(
                    boost::make_tuple(cm->first, FieldEntity::getHiBitNumberUId(
                                                     BITFIELDID_SIZE - 1)));
 
            auto sub_inv_mat =
                matrix_range(invMat, get_range(rbi), get_range(cbi));
            bM.resize(sub_inv_mat.size1(), sub_inv_mat.size2());
            noalias(bM) = sub_inv_mat;
 
            for (auto a_lo = a_lo_tmp; a_lo != a_hi; ++a_lo) {
#ifndef NDEBUG
              if (block_indexing.find((*a_lo)->uidRow) != block_indexing.end())
                SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                        "Wrong a_lo->uidRow");
#endif
 
              auto &a = (*a_lo)->getMat();
              abM.resize(a.size1(), bM.size2(), false);
              abM.clear();
              cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
 
                          a.size1(), bM.size2(), a.size2(), 1.,
 
                          &*a.data().begin(), a.size2(),
 
                          &*bM.data().begin(), bM.size2(), 0.,
 
                          &*abM.data().begin(), abM.size2());
 
              for (auto c_lo = c_lo_tmp; c_lo != c_hi; ++c_lo) {
#ifndef NDEBUG
                if (block_indexing.find((*c_lo)->uidCol) !=
                    block_indexing.end())
                  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                          "Wrong a_lo->uidRow");
#endif
 
                auto &c = (*c_lo)->getMat();
 
                abcM.resize(abM.size1(), c.size2(), false);
                abcM.clear();
 
                cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                            abM.size1(), c.size2(), abM.size2(), -1.,
 
                            &*abM.data().begin(), abM.size2(),
 
                            &*c.data().begin(), c.size2(), 0.,
 
                            &*abcM.data().begin(), abcM.size2());
 
                CHKERR assemble_schur(abcM, (*a_lo)->uidRow, (*c_lo)->uidCol,
                                      &((*a_lo)->getRowInd()),
                                      &((*c_lo)->getColInd()));
 
                if (symSchur && rp_uid_it != cp_uid_it) {
                  abcM = trans(abcM);
                  CHKERR assemble_schur(abcM, (*c_lo)->uidCol, (*a_lo)->uidRow,
                                        &((*c_lo)->getColInd()),
                                        &((*a_lo)->getRowInd()));
                }
              }
            }
          }
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  // Assemble Schur complements
  CHKERR assemble_dense_blocks();
 
#ifndef NDEBUG
  if constexpr (debug_schur)
    MOFEM_LOG("SELF", Sev::noisy) << "Schur assemble done";
#endif
 
#ifndef NDEBUG
  PetscLogEventEnd(SchurEvents::MOFEM_EVENT_opSchurAssembleEnd, 0, 0, 0, 0);
#endif
 
  MoFEMFunctionReturn(0);
}