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. More...

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 74 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 76 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 510 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 521 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_field_name = [&](auto uid) {
     return OP::getPtrFE()->mField.get_field_name(field_bit_from_bit_number(
         FieldEntity::getFieldBitNumberFromUniqueId(uid)));
   };
  
   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 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 [block_list, block_indexing, block_mat_size] =
         get_block_indexing(a00_uids);
  
     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) {
       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);
 }