49 return ::MatSetValues(
51 S, row_ind.size(), &*row_ind.begin(), col_ind.size(), &*col_ind.begin(),
52 &*
m.data().begin(), iora
72 boost::shared_ptr<std::vector<unsigned char>> marker_ptr,
79 boost::shared_ptr<std::vector<unsigned char>>
markerPtr;
87template <
typename OP_SCHUR_ASSEMBLE_BASE>
105 std::vector<std::string> fields_name,
106 std::vector<boost::shared_ptr<Range>> field_ents,
110 bool sym_schur,
bool symm_op);
113 template <
typename I>
161struct SchurElemMats :
public boost::enable_shared_from_this<SchurElemMats> {
187 template <
typename OP_SCHUR_ASSEMBLE_BASE>
202 member<SchurElemMats, const UId, &SchurElemMats::uidRow>,
203 member<SchurElemMats, const UId, &SchurElemMats::uidCol>
212 static boost::ptr_vector<MatrixDouble>
locMats;
279 const_mem_fun<Indexes, UId, &Indexes::getRowUId>,
280 const_mem_fun<Indexes, UId, &Indexes::getColUId>>
286 const_mem_fun<Indexes, UId, &Indexes::getFEUId>
292 const_mem_fun<Indexes, UId, &Indexes::getRowUId>
298 const_mem_fun<Indexes, UId, &Indexes::getColUId>
344 : iDX(idx), uidRow(uid_row), uidCol(uid_col) {}
358 auto get_row_indices = [&]() ->
const VectorInt & {
360 if (
auto stored_data_ptr =
362 return stored_data_ptr->entityIndices;
368 const auto &row_ind = get_row_indices();
371 const auto nb_rows = row_ind.size();
372 const auto nb_cols = col_ind.size();
375 if (mat.size1() != nb_rows) {
377 "Wrong mat size %ld != %ld", mat.size1(), nb_rows);
379 if (mat.size2() != nb_cols) {
381 "Wrong mat size %ld != %ld", mat.size2(), nb_cols);
386 auto get_uid = [](
auto &data) {
387 if (data.getFieldEntities().size() == 1) {
389 return data.getFieldEntities()[0]->getLocalUniqueId();
398 auto &uid0 = data.getFieldEntities()[0]->getLocalUniqueId();
404 for (
auto i = 1;
i < data.getFieldEntities().size(); ++
i) {
411 data.getFieldEntities()[
i]->getLocalUniqueId())
425 auto uid_row = get_uid(row_data);
426 auto uid_col = get_uid(col_data);
430 .find(boost::make_tuple(uid_row, uid_col));
462 auto asmb = [&](
auto &sm) {
463 sm.resize(nb_rows, nb_cols,
false);
467 asmb((*p.first)->getMat());
469 auto add_indices = [](
auto &storage,
auto &ind) {
470 storage.resize(ind.size(),
false);
471 noalias(storage) = ind;
474 add_indices((*p.first)->getRowInd(), row_ind);
475 add_indices((*p.first)->getColInd(), col_ind);
480 auto asmb = [&](
auto &sm) {
484 if (sm.size1() != nb_rows) {
486 "Wrong mat or storage size %ld != %ld", sm.size1(), nb_rows);
488 if (sm.size2() != nb_cols) {
490 "Wrong mat or storage size %ld != %ld", sm.size2(), nb_cols);
503 "Assembly type not implemented");
508 CHKERR asmb((*it)->getMat());
528 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble begin";
537 boost::shared_ptr<std::vector<unsigned char>> marker_ptr,
double diag_val)
547 auto prb_ptr = fe_ptr->problemPtr;
550 auto find_loc_dof_index = [&](
auto glob) {
551 auto it = dof.find(glob);
552 if (it != dof.end()) {
553 return (*it)->getPetscLocalDofIdx();
555 MOFEM_LOG(
"SELF", Sev::error) <<
"Global dof " << glob;
558 "Global dof index not found in local numered dofs storage");
564 auto zero_ind = [&](
auto &ind) {
565 for (
auto i = 0;
i < ind.size(); ++
i) {
567 if ((*
markerPtr)[find_loc_dof_index(ind[
i])]) {
574 zero_ind(s->getRowInd());
575 zero_ind(s->getColInd());
582template <
typename OP_SCHUR_ASSEMBLE_BASE>
584 std::vector<std::string> fields_name,
585 std::vector<boost::shared_ptr<Range>> field_ents,
588 :
OP(
NOSPACE,
OP::OPSPACE, symm_op), fieldsName(fields_name),
589 fieldEnts(field_ents), schurAO(schur_ao), schurMat(schur_mat),
590 symSchur(sym_schur) {}
592template <
typename OP_SCHUR_ASSEMBLE_BASE>
605 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble begin -> end";
608 auto get_a00_uids = [&]() {
609 auto get_field_bit = [&](
auto &name) {
610 return OP::getPtrFE()->mField.get_field_bit_number(name);
613 std::vector<std::pair<UId, UId>> a00_uids;
614 a00_uids.reserve(fieldsName.size());
615 for (
auto ss = 0; ss != fieldsName.size(); ++ss) {
616 auto field_bit = get_field_bit(fieldsName[ss]);
617 auto row_ents = fieldEnts[ss];
619 for (
auto p = row_ents->pair_begin(); p != row_ents->pair_end(); ++p) {
624 a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
631 a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
638 auto get_field_name = [&](
auto uid) {
642 auto list_storage = [&](
auto &storage) {
645 for (
auto &p : storage) {
647 <<
"List schur storage: " <<
i <<
" " << p->iDX <<
": "
648 << get_field_name(p->uidRow) <<
" " << get_field_name(p->uidCol)
649 <<
" : " << p->getMat().size1() <<
" " << p->getMat().size2();
656 auto assemble_dense_blocks = [&]() {
657 using matrix_range = ublas::matrix_range<MatrixDouble>;
658 using range = ublas::range;
662 auto assemble_schur = [
this](
auto &
m,
auto &uid_row,
auto &uid_col,
663 auto *row_ind_ptr,
auto *col_ind_ptr) {
668 if (
auto ierr = this->assembleSchurMat(
670 schurMat, uid_row, *row_ind_ptr, uid_col, *col_ind_ptr,
m,
675 auto field_ents = OP::getPtrFE()->mField.get_field_ents();
676 auto row_ent_it = field_ents->find(uid_row);
677 auto col_ent_it = field_ents->find(uid_col);
679 if (row_ent_it != field_ents->end())
681 <<
"Assemble row entity: " << (*row_ent_it)->getName() <<
" "
682 << (*col_ent_it)->getEntTypeName() <<
" side "
683 << (*row_ent_it)->getSideNumber();
684 if (col_ent_it != field_ents->end())
686 <<
"Assemble col entity: " << (*col_ent_it)->getName() <<
" "
687 << (*col_ent_it)->getEntTypeName() <<
" side "
688 << (*col_ent_it)->getSideNumber();
697 auto a00_uids = get_a00_uids();
699 auto get_block_indexing = [&](
auto &a00_uids) {
701 std::vector<const SchurElemMats *> block_list;
702 block_list.reserve(storage.size());
703 for (
auto &rp_uid : a00_uids) {
704 auto [rlo_uid, rhi_uid] = rp_uid;
705 for (
auto &cp_uid : a00_uids) {
706 auto [clo_uid, chi_uid] = cp_uid;
709 storage.template get<SchurElemMats::uid_mi_tag>().lower_bound(
710 boost::make_tuple(rlo_uid, clo_uid));
712 storage.template get<SchurElemMats::uid_mi_tag>().upper_bound(
713 boost::make_tuple(rhi_uid, chi_uid));
715 for (; it != hi_it; ++it) {
716 if ((*it)->uidRow >= rlo_uid && (*it)->uidRow < rhi_uid &&
717 (*it)->uidCol >= clo_uid && (*it)->uidCol < chi_uid) {
718 block_list.push_back(*it);
725 std::map<UId, std::pair<size_t, const VectorInt *>>
727 for (
auto d : block_list) {
728 if (block_indexing.find(d->uidRow) == block_indexing.end()) {
729 block_indexing[d->uidRow] =
730 std::make_pair(d->getRowInd().size(), &(d->getRowInd()));
735 int mat_block_size = 0;
736 for (
auto &p_uid : a00_uids) {
737 auto [lo_uid, hi_uid] = p_uid;
738 auto lo = block_indexing.lower_bound(lo_uid);
739 auto up = block_indexing.upper_bound(hi_uid);
740 for (; lo != up; ++lo) {
741 lo->second.first = mat_block_size;
742 mat_block_size += lo->second.second->size();
746 return std::make_tuple(block_list, block_indexing, mat_block_size);
749 auto get_schur_block_list = [&](
auto &block_indexing) {
750 std::vector<const SchurElemMats *> block_list;
751 block_list.reserve(storage.size());
752 for (
auto &s : storage) {
753 if (block_indexing.find(s->uidRow) == block_indexing.end() &&
754 block_indexing.find(s->uidCol) == block_indexing.end()) {
755 block_list.push_back(s);
761 auto bi = get_block_indexing(a00_uids);
762 auto &[block_list, block_indexing, block_mat_size] = bi;
764 if (block_mat_size == 0) {
767 for (
auto &s : storage) {
768 CHKERR AOApplicationToPetsc(schurAO, s->getRowInd().size(),
769 &*s->getRowInd().begin());
770 CHKERR AOApplicationToPetsc(schurAO, s->getColInd().size(),
771 &*s->getColInd().begin());
775 for (
auto &s : storage) {
776 auto &
m = s->getMat();
777 CHKERR assemble_schur(
m, s->uidRow, s->uidCol, &(s->getRowInd()),
783 blockMat.resize(block_mat_size, block_mat_size,
false);
786 auto get_range = [](
auto &bi) {
787 return range(bi.first, bi.first + bi.second->size());
790 for (
auto &s : block_list) {
791 auto &
m = s->getMat();
793 if (block_indexing.find(s->uidRow) == block_indexing.end())
795 if (block_indexing.find(s->uidCol) == block_indexing.end())
799 auto &rbi = block_indexing.at(s->uidRow);
800 auto &cbi = block_indexing.at(s->uidCol);
802 auto sub_mat = matrix_range(blockMat, get_range(rbi), get_range(cbi));
806 auto get_zeroed_indices = [&](
auto extractor_uid,
auto extractor_ind) {
807 auto &[block_list, block_indexing, block_mat_size] = bi;
808 std::vector<int> zeroed_indices;
809 zeroed_indices.reserve(block_mat_size);
810 for (
auto &s : block_list) {
811 auto &bi = block_indexing.at(extractor_uid(s));
812 auto &ind = extractor_ind(s);
813 for (
auto it = ind.begin(); it != ind.end(); ++it) {
815 auto idx = bi.first + std::distance(ind.begin(), it);
816 zeroed_indices.push_back(idx);
820 std::sort(zeroed_indices.begin(), zeroed_indices.end());
821 auto it = std::unique(zeroed_indices.begin(), zeroed_indices.end());
822 zeroed_indices.resize(std::distance(zeroed_indices.begin(), it));
823 return zeroed_indices;
825 auto zero_rows = get_zeroed_indices(
826 [](
auto &s) {
return s->uidRow; },
827 [](
auto &s) ->
VectorInt & {
return s->getRowInd(); });
828 auto zero_cols = get_zeroed_indices(
829 [](
auto &s) {
return s->uidCol; },
830 [](
auto &s) ->
VectorInt & {
return s->getColInd(); });
832 for (
auto i : zero_rows) {
833 for (
auto j = 0;
j != blockMat.size2(); ++
j) {
837 for (
auto j : zero_cols) {
838 for (
auto i = 0;
i != blockMat.size1(); ++
i) {
842 for (
auto i : zero_rows) {
846 CHKERR I::invertMat(blockMat, invMat);
849 for (
auto &s : block_list) {
850 auto it = storage.template get<SchurElemMats::uid_mi_tag>().find(
851 boost::make_tuple(s->uidRow, s->uidCol));
852 storage.template get<SchurElemMats::uid_mi_tag>().erase(it);
857 for (
auto &s : storage) {
858 if (block_indexing.find(s->uidRow) == block_indexing.end()) {
859 CHKERR AOApplicationToPetsc(schurAO, s->getRowInd().size(),
860 &*s->getRowInd().begin());
862 if (block_indexing.find(s->uidCol) == block_indexing.end()) {
863 CHKERR AOApplicationToPetsc(schurAO, s->getColInd().size(),
864 &*s->getColInd().begin());
869 auto schur_block_list = get_schur_block_list(block_indexing);
870 for (
auto &s : schur_block_list) {
871 auto &
m = s->getMat();
872 CHKERR assemble_schur(
m, s->uidRow, s->uidCol, &(s->getRowInd()),
875 for (
auto &s : schur_block_list) {
876 auto it = storage.template get<SchurElemMats::uid_mi_tag>().find(
877 boost::make_tuple(s->uidRow, s->uidCol));
878 storage.template get<SchurElemMats::uid_mi_tag>().erase(it);
880 schur_block_list.clear();
884 auto rp_uid_it = a00_uids.begin(); rp_uid_it != a00_uids.end();
888 auto [rlo_uid, rhi_uid] = *rp_uid_it;
889 for (
auto rm = block_indexing.lower_bound(rlo_uid);
890 rm != block_indexing.upper_bound(rhi_uid); ++rm) {
891 auto &rbi = rm->second;
894 storage.template get<SchurElemMats::col_mi_tag>().lower_bound(
897 storage.template get<SchurElemMats::col_mi_tag>().upper_bound(
902 auto cp_uid_it = (symSchur) ? rp_uid_it : a00_uids.begin();
903 cp_uid_it != a00_uids.end(); ++cp_uid_it
906 auto [clo_uid, chi_uid] = *cp_uid_it;
907 for (
auto cm = block_indexing.lower_bound(clo_uid);
908 cm != block_indexing.upper_bound(chi_uid); ++cm) {
909 auto &cbi = cm->second;
912 storage.template get<SchurElemMats::uid_mi_tag>().lower_bound(
913 boost::make_tuple(cm->first, 0));
915 storage.template get<SchurElemMats::uid_mi_tag>().upper_bound(
920 matrix_range(invMat, get_range(rbi), get_range(cbi));
921 bM.resize(sub_inv_mat.size1(), sub_inv_mat.size2());
922 noalias(bM) = sub_inv_mat;
924 for (
auto a_lo = a_lo_tmp; a_lo != a_hi; ++a_lo) {
926 if (block_indexing.find((*a_lo)->uidRow) != block_indexing.end())
928 "Wrong a_lo->uidRow");
931 auto &
a = (*a_lo)->getMat();
932 abM.resize(
a.size1(), bM.size2(),
false);
934 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
936 a.size1(), bM.size2(),
a.size2(), 1.,
938 &*
a.data().begin(),
a.size2(),
940 &*bM.data().begin(), bM.size2(), 0.,
942 &*abM.data().begin(), abM.size2());
944 for (
auto c_lo = c_lo_tmp; c_lo != c_hi; ++c_lo) {
946 if (block_indexing.find((*c_lo)->uidCol) !=
947 block_indexing.end())
949 "Wrong a_lo->uidRow");
952 auto &
c = (*c_lo)->getMat();
954 abcM.resize(abM.size1(),
c.size2(),
false);
957 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
958 abM.size1(),
c.size2(), abM.size2(), -1.,
960 &*abM.data().begin(), abM.size2(),
962 &*
c.data().begin(),
c.size2(), 0.,
964 &*abcM.data().begin(), abcM.size2());
966 CHKERR assemble_schur(abcM, (*a_lo)->uidRow, (*c_lo)->uidCol,
967 &((*a_lo)->getRowInd()),
968 &((*c_lo)->getColInd()));
970 if (symSchur && rp_uid_it != cp_uid_it) {
972 CHKERR assemble_schur(abcM, (*c_lo)->uidCol, (*a_lo)->uidRow,
973 &((*c_lo)->getColInd()),
974 &((*a_lo)->getRowInd()));
987 CHKERR assemble_dense_blocks();
991 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble done";
1006 const auto nb =
m.size1();
1008 if (nb !=
m.size2()) {
1010 "It should be square matrix %ld != %ld", nb,
m.size2());
1014 inv.resize(nb, nb,
false);
1016 m.resize(2 * nb, 2 * nb,
false);
1023 info =
lapack_dgetrf(nb, nb, &*inv.data().begin(), nb, &*ipiv.begin());
1026 "lapack error info = %d", info);
1027 info =
lapack_dgetri(nb, &*inv.data().begin(), nb, &*ipiv.begin(),
1028 &*
m.data().begin(),
m.data().size());
1031 "lapack error info = %d", info);
1042 const auto nb =
m.size1();
1044 if (nb !=
m.size2()) {
1046 "It should be square matrix %ld != %ld", nb,
m.size2());
1050 inv.resize(nb, nb,
false);
1056 info =
lapack_dgetrf(nb, nb, &*inv.data().begin(), nb, &*ipiv.begin());
1059 "lapack error info = %d", info);
1060 info =
lapack_dgetri(nb, &*inv.data().begin(), nb, &*ipiv.begin(),
1061 &*
m.data().begin(),
m.data().size());
1064 "lapack error info = %d", info);
1073 return doWorkImpl<SchurDSYSV>(side, type, data);
1079 return doWorkImpl<SchurDGESV>(side, type, data);
1089 auto cmp_uid_lo = [](
const boost::weak_ptr<FieldEntity> &
a,
const UId &b) {
1090 if (
auto a_ptr =
a.lock()) {
1091 if (a_ptr->getLocalUniqueId() < b)
1100 auto cmp_uid_hi = [](
const UId &b,
const boost::weak_ptr<FieldEntity> &
a) {
1101 if (
auto a_ptr =
a.lock()) {
1102 if (b < a_ptr->getLocalUniqueId())
1112 auto get_uid_pair = [](
const auto &field_id) {
1117 return std::make_pair(lo_uid, hi_uid);
1122 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(),
1124 auto hi = std::upper_bound(field_ents.begin(), field_ents.end(),
1126 return std::make_pair(lo, hi);
1131 auto row_extractor = [](
auto &e) {
return e->entityCacheRowDofs; };
1132 auto col_extractor = [](
auto &e) {
return e->entityCacheColDofs; };
1134 auto extract_data = [](
auto &&its,
auto extractor) {
1135 std::vector<std::tuple<UId, int, int, int>> data;
1136 data.reserve(std::distance(its.first, its.second));
1138 for (; its.first != its.second; ++its.first) {
1139 if (
auto e = its.first->lock()) {
1140 if (
auto cache = extractor(e).lock()) {
1142 for (
auto dof = cache->loHi[0]; dof != cache->loHi[1]; ++dof) {
1143 if ((*dof)->getPetscGlobalDofIdx() != -1)
1146 auto uid = e->getLocalUniqueId();
1149 auto glob = (*cache->loHi[0])->getPetscGlobalDofIdx();
1150 auto loc = (*cache->loHi[0])->getPetscLocalDofIdx();
1151 while (glob == -1 && cache->loHi[0] != cache->loHi[1]) {
1153 glob = (*cache->loHi[0])->getPetscGlobalDofIdx();
1154 loc = (*cache->loHi[0])->getPetscLocalDofIdx();
1156 data.emplace_back(uid, glob, nb_dofs, loc);
1160 for (
auto lo = cache->loHi[0]; lo != cache->loHi[1]; ++lo) {
1161 auto glob = (*lo)->getPetscGlobalDofIdx();
1164 "Wrong global index");
1170 data.emplace_back(uid, -1, 0, -1);
1178 auto data_ptr = boost::make_shared<BlockStructure>();
1182 auto fe_method = boost::shared_ptr<MoFEM::FEMethod>(
new MoFEM::FEMethod());
1184 for (
auto &d : schur_fe_op_vec) {
1187 auto get_bit_numbers = [&d](
auto op) {
1188 std::vector<FieldBitNumber> bit_numbers(d.second.size());
1189 std::transform(d.second.begin(), d.second.end(), bit_numbers.begin(), op);
1194 auto row_bit_numbers = get_bit_numbers(
1195 [&](
auto &p) {
return m_field_ptr->get_field_bit_number(p.first); });
1197 auto col_bit_numbers = get_bit_numbers(
1198 [&](
auto &p) {
return m_field_ptr->get_field_bit_number(p.second); });
1200 fe_method->preProcessHook = []() {
return 0; };
1201 fe_method->postProcessHook = []() {
return 0; };
1202 fe_method->operatorHook = [&]() {
1205 auto fe_uid = fe_method->numeredEntFiniteElementPtr->getLocalUniqueId();
1207 for (
auto f = 0; f != row_bit_numbers.size(); ++f) {
1210 extract_data(get_it_pair(fe_method->getRowFieldEnts(),
1211 get_uid_pair(row_bit_numbers[f])),
1214 extract_data(get_it_pair(fe_method->getColFieldEnts(),
1215 get_uid_pair(col_bit_numbers[f])),
1218 for (
auto &r : row_data) {
1219 auto [r_uid, r_glob, r_nb_dofs, r_loc] = r;
1220 for (
auto &
c : col_data) {
1221 auto [c_uid, c_glob, c_nb_dofs, c_loc] =
c;
1223 r_uid, c_uid, fe_uid, r_glob, c_glob, r_nb_dofs, c_nb_dofs,
1233 m_field_ptr->get_comm_size());
1238 for (
auto &
v : data_ptr->blockIndex.get<0>()) {
1239 v.getMatShift() = mem_size;
1240 mem_size +=
v.getNbCols() *
v.getNbRows();
1243 std::vector<double> tmp;
1244 if (mem_size > tmp.max_size())
1247 data_ptr->dataBlocksPtr =
1248 boost::make_shared<std::vector<double>>(mem_size, 0);
1252 CHKERR VecSetDM(ghost_x, PETSC_NULLPTR);
1253 CHKERR VecSetDM(ghost_y, PETSC_NULLPTR);
1255 data_ptr->ghostX = ghost_x;
1256 data_ptr->ghostY = ghost_y;
1262 Mat mat, Vec x, Vec y, InsertMode iora,
1264 int(DiagBlockIndex::BlockIndex::nth_index<0>::type::iterator)>
1266 boost::shared_ptr<std::vector<double>> data_blocks_ptr,
1267 bool multiply_by_preconditioner);
1272static PetscErrorCode
mult(Mat mat, Vec x, Vec y) {
1274 CHKERR MatShellGetContext(mat, (
void **)&ctx);
1276 mat, x, y, INSERT_VALUES,
1278 [](DiagBlockIndex::BlockIndex::nth_index<0>::type::iterator it) {
1279 return it->getMatShift();
1286 CHKERR MatShellGetContext(mat, (
void **)&ctx);
1288 mat, x, y, ADD_VALUES,
1290 [](DiagBlockIndex::BlockIndex::nth_index<0>::type::iterator it) {
1291 return it->getMatShift();
1296static PetscErrorCode
solve(Mat mat, Vec x, Vec y) {
1304 const PetscInt rows[], PetscScalar diag,
1309 CHKERR MatShellGetContext(A, (
void **)&ctx);
1313 using BlockIndexView = multi_index_container<
1335 BlockIndexView view;
1336 auto hint = view.get<0>().end();
1338 hint = view.insert(hint, &
v);
1341 const int *ptr = &rows[0];
1346 CHKERR PetscObjectGetComm((PetscObject)A, &comm);
1348 MPI_Comm_size(comm, &size);
1354 CHKERR ISGetSize(is_local, &is_nb_rows);
1357 CHKERR ISView(is_local, PETSC_VIEWER_STDOUT_WORLD);
1360 CHKERR ISGetIndices(is_local, &ptr);
1364 CHKERR MatGetLocalSize(A, &loc_m, &loc_n);
1366 for (
auto n = 0;
n != is_nb_rows; ++
n) {
1368 auto rlo = view.get<0>().lower_bound(row);
1369 auto rhi = view.get<0>().end();
1370 for (; rlo != rhi; ++rlo) {
1371 auto r_shift = row - (*rlo)->getRow();
1372 if (r_shift >= 0 && r_shift < (*rlo)->getNbRows()) {
1374 for (
auto i = 0;
i != (*rlo)->getNbCols(); ++
i) {
1375 ptr[
i + r_shift * (*rlo)->getNbCols()] = 0;
1377 }
else if ((*rlo)->getRow() + (*rlo)->getNbRows() > row) {
1383 for (
auto n = 0;
n != is_nb_rows; ++
n) {
1385 auto clo = view.get<1>().lower_bound(col);
1386 auto chi = view.get<1>().end();
1387 for (; clo != chi; ++clo) {
1388 auto c_shift = col - (*clo)->getCol();
1389 if (c_shift >= 0 && c_shift < (*clo)->getNbCols()) {
1392 for (
auto i = 0;
i != (*clo)->getNbRows(); ++
i) {
1393 ptr[c_shift +
i * (*clo)->getNbCols()] = 0;
1399 (*clo)->getRow() == (*clo)->getCol() &&
1400 (*clo)->getLocRow() < loc_m && (*clo)->getLocCol() < loc_n
1403 auto r_shift = col - (*clo)->getCol();
1404 if (r_shift >= 0 && r_shift < (*clo)->getNbRows()) {
1405 ptr[c_shift + r_shift * (*clo)->getNbCols()] = diag;
1408 }
else if ((*clo)->getCol() + (*clo)->getNbCols() > col) {
1415 CHKERR ISRestoreIndices(is_local, &ptr);
1426 CHKERR MatShellGetContext(
m, (
void **)&ctx);
1437 CHKERR MatShellSetManageScalingShifts(mat_raw);
1438 CHKERR MatShellSetOperation(mat_raw, MATOP_MULT, (
void (*)(
void))
mult);
1439 CHKERR MatShellSetOperation(mat_raw, MATOP_MULT_ADD,
1441 CHKERR MatShellSetOperation(mat_raw, MATOP_SOLVE, (
void (*)(
void))
solve);
1442 CHKERR MatShellSetOperation(mat_raw, MATOP_SOLVE_ADD,
1444 CHKERR MatShellSetOperation(mat_raw, MATOP_ZERO_ENTRIES,
1446 CHKERR MatShellSetOperation(mat_raw, MATOP_ZERO_ROWS_COLUMNS,
1453 boost::shared_ptr<BlockStructure> data) {
1456 auto nb_local = problem_ptr->nbLocDofsRow;
1457 auto nb_global = problem_ptr->nbDofsRow;
1460 if (nb_local != problem_ptr->nbLocDofsCol) {
1462 <<
"Wrong size " << nb_local <<
" != " << problem_ptr->nbLocDofsCol;
1464 "nb. cols is inconsistent with nb. rows");
1466 if (nb_global != problem_ptr->nbDofsCol) {
1468 <<
"Wrong size " << nb_global <<
" != " << problem_ptr->nbDofsCol;
1470 "nb. cols is inconsistent with nb. rows");
1475 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
1478 CHKERR MatCreateShell(comm, nb_local, nb_local, nb_global, nb_global,
1479 (
void *)data.get(), &mat_raw);
1489 Mat mat, Vec x, Vec y, InsertMode iora,
1492 int(DiagBlockIndex::BlockIndex::nth_index<0>::type::iterator)>
1495 boost::shared_ptr<std::vector<double>> data_blocks_ptr,
1496 bool multiply_by_preconditioner) {
1499 CHKERR MatShellGetContext(mat, (
void **)&ctx);
1504 CHKERR VecGetLocalSize(x, &x_loc_size);
1506 CHKERR VecGetLocalSize(y, &y_loc_size);
1508 Vec ghost_x = ctx->
ghostX;
1509 Vec ghost_y = ctx->
ghostY;
1511 CHKERR VecCopy(x, ghost_x);
1515 CHKERR VecGhostGetLocalForm(ghost_y, &loc_ghost_y);
1517 CHKERR VecGetLocalSize(loc_ghost_y, &nb_y);
1518 CHKERR VecGetArray(loc_ghost_y, &y_array);
1519 for (
auto i = 0;
i != nb_y; ++
i)
1521 CHKERR VecRestoreArray(loc_ghost_y, &y_array);
1522 CHKERR VecGhostRestoreLocalForm(ghost_y, &loc_ghost_y);
1524 auto mult = [&](
int low_x,
int hi_x,
int low_y,
int hi_y) {
1527 const double *x_array;
1529 CHKERR VecGhostGetLocalForm(ghost_x, &loc_ghost_x);
1530 CHKERR VecGetArrayRead(loc_ghost_x, &x_array);
1534 CHKERR VecGhostGetLocalForm(ghost_y, &loc_ghost_y);
1536 CHKERR VecGetLocalSize(loc_ghost_y, &nb_y);
1537 CHKERR VecGetArray(loc_ghost_y, &y_array);
1539 double *block_ptr = &*data_blocks_ptr->begin();
1544 if (it->getLocRow() < low_y || it->getLocRow() >= hi_y ||
1545 it->getLocCol() < low_x || it->getLocCol() >= hi_x) {
1550 auto nb_rows = it->getNbRows();
1551 auto nb_cols = it->getNbCols();
1552 auto x_ptr = &x_array[it->getLocCol()];
1553 auto y_ptr = &y_array[it->getLocRow()];
1554 auto ptr = &block_ptr[shift_extractor(it)];
1557 cblas_dgemv(CblasRowMajor, CblasNoTrans, nb_rows, nb_cols, 1.0, ptr,
1558 nb_cols, x_ptr, 1, 1.0, y_ptr, 1);
1560 for (
auto r = 0; r != nb_rows; ++r) {
1561 for (
auto c = 0;
c != nb_cols; ++
c) {
1562 y_ptr[r] += ptr[r * nb_cols +
c] * x_ptr[
c];
1573 "No parentBlockStructurePtr");
1581 if (it->getLocRow() < low_y || it->getLocRow() >= hi_y ||
1582 it->getLocCol() < low_x || it->getLocCol() >= hi_x) {
1587 if (it->getMatShift() != -1) {
1588 auto nb_rows = it->getNbRows();
1589 auto nb_cols = it->getNbCols();
1590 auto x_ptr = &x_array[it->getLocCol()];
1591 auto y_ptr = &y_array[it->getLocRow()];
1592 auto ptr = &preconditioner_ptr[it->getMatShift()];
1594 cblas_dgemv(CblasRowMajor, CblasNoTrans, nb_rows, nb_cols, 1.0, ptr,
1595 nb_cols, x_ptr, 1, 1.0, y_ptr, 1);
1597 for (
auto r = 0; r != nb_rows; ++r) {
1598 for (
auto c = 0;
c != nb_cols; ++
c) {
1599 y_ptr[r] += ptr[r * nb_cols +
c] * x_ptr[
c];
1609 CHKERR VecRestoreArrayRead(loc_ghost_x, &x_array);
1610 CHKERR VecRestoreArray(loc_ghost_y, &y_array);
1611 CHKERR VecGhostRestoreLocalForm(ghost_x, &loc_ghost_x);
1612 CHKERR VecGhostRestoreLocalForm(ghost_y, &loc_ghost_y);
1616 constexpr auto max_int = std::numeric_limits<int>::max();
1617 CHKERR VecGhostUpdateBegin(ghost_x, INSERT_VALUES, SCATTER_FORWARD);
1619 CHKERR VecGhostUpdateEnd(ghost_x, INSERT_VALUES, SCATTER_FORWARD);
1620 CHKERR mult(x_loc_size, max_int, 0, max_int);
1622 CHKERR VecGhostUpdateBegin(ghost_y, ADD_VALUES, SCATTER_REVERSE);
1623 CHKERR VecGhostUpdateEnd(ghost_y, ADD_VALUES, SCATTER_REVERSE);
1627 CHKERR VecCopy(ghost_y, y);
1630 CHKERR VecAXPY(y, 1., ghost_y);
1638 auto print_norm = [&](
auto msg,
auto y) {
1641 CHKERR VecNorm(y, NORM_2, &norm);
1643 CHKERR VecGetLocalSize(y, &nb_loc_y);
1645 CHKERR VecGetSize(y, &nb_y);
1647 << msg <<
" " << nb_y <<
" " << nb_loc_y <<
" norm " << norm;
1653 print_norm(
"mult_schur_block_shell insert x", x);
1654 print_norm(
"mult_schur_block_shell insert y", y);
1657 print_norm(
"mult_schur_block_shell add x", x);
1658 print_norm(
"mult_schur_block_shell add y", y);
1674 using matrix_range = ublas::matrix_range<MatrixDouble>;
1675 using range = ublas::range;
1678 CHKERR MatShellGetContext(mat, (
void **)&ctx);
1685 "No preconditionerBlocksPtr");
1688 if (iora == INSERT_VALUES) {
1689 CHKERR VecZeroEntries(x);
1693 CHKERR VecGetArray(x, &x_array);
1694 const double *y_array;
1695 CHKERR VecGetArrayRead(y, &y_array);
1697 std::vector<const DiagBlockIndex::Indexes *> blocks;
1702 auto &block_index = ctx->
blockIndex.get<1>();
1703 for (
auto it = block_index.begin(); it != block_index.end();) {
1707 auto last_uid = it->getFEUId();
1708 while (it != block_index.end() && it->getFEUId() == last_uid) {
1709 blocks.push_back(&*it);
1714 std::map<UId, std::tuple<int, int, int>> block_indexing;
1715 for (
auto b : blocks) {
1716 if (block_indexing.find(b->getRowUId()) == block_indexing.end()) {
1717 block_indexing[b->getRowUId()] =
1718 std::make_tuple(b->getNbRows(), b->getNbRows(), b->getLocRow());
1721 std::sort(blocks.begin(), blocks.end(), [](
auto a,
auto b) {
1722 if (a->getRowUId() == b->getRowUId())
1723 return a->getColUId() < b->getColUId();
1725 return a->getRowUId() < b->getRowUId();
1729 int mat_block_size = 0;
1730 for (
auto &b : block_indexing) {
1731 auto &[index, size, loc] = b.second;
1732 index = mat_block_size;
1733 mat_block_size += size;
1736 std::vector<std::tuple<int, int, int, int, int>> block_data;
1737 block_data.reserve(blocks.size());
1738 for (
auto s : blocks) {
1739 auto ruid = s->getRowUId();
1740 auto cuid = s->getColUId();
1741 auto &rbi = block_indexing.at(ruid);
1742 auto &cbi = block_indexing.at(cuid);
1743 if (
auto shift = s->getMatShift(); shift != -1) {
1744 block_data.push_back(std::make_tuple(
1748 s->getNbRows(), s->getNbCols(),
1750 get<0>(rbi), get<0>(cbi))
1755 block_mat.resize(mat_block_size, mat_block_size,
false);
1757 for (
auto &bd : block_data) {
1758 auto &[shift, nb_rows, nb_cols, ridx, cidx] = bd;
1760 for (
auto i = 0;
i != nb_rows; ++
i, ptr += nb_cols) {
1762 cblas_dcopy(nb_cols, ptr, 1, sub_ptr, 1);
1766 for (
auto &bd : block_data) {
1767 auto &[shift, nb_rows, nb_cols, ridx, cidx] = bd;
1769 for (
auto i = 0;
i != nb_rows; ++
i, ptr += nb_cols) {
1771 cblas_daxpy(nb_cols, 1., ptr, 1, sub_ptr, 1);
1777 block_y.resize(mat_block_size,
false);
1780 for (
auto &b : block_indexing) {
1781 auto [index, size, loc] = b.second;
1782 auto ptr = &y_array[loc];
1783 cblas_dcopy(size, ptr, 1, &block_y[index], 1);
1787 constexpr int nrhs = 1;
1788 ipiv.resize(mat_block_size,
false);
1790 mat_block_size, &*ipiv.data().begin(),
1791 &*block_y.data().begin(), mat_block_size);
1794 <<
"error lapack solve dgesv info = " << info;
1797 "error lapack solve dgesv info = %d", info);
1800 for (
auto &b : block_indexing) {
1801 auto [index, size, loc] = b.second;
1802 auto ptr = &x_array[loc];
1803 cblas_daxpy(size, 1.0, &block_y[index], 1, ptr, 1);
1807 CHKERR VecRestoreArray(x, &x_array);
1808 CHKERR VecRestoreArrayRead(y, &y_array);
1818 std::vector<std::string> fields_name,
1819 std::vector<boost::shared_ptr<Range>> field_ents,
1821 using matrix_range = ublas::matrix_range<MatrixDouble>;
1822 using range = ublas::range;
1825 CHKERR MatShellGetContext(
B, (
void **)&ctx);
1827 constexpr bool debug =
false;
1838 "No preconditionerBlocksPtr");
1841 std::vector<std::pair<UId, UId>> a00_fields_uids;
1842 a00_fields_uids.reserve(fields_name.size());
1843 auto it_field_ents = field_ents.begin();
1844 if (fields_name.size() != field_ents.size()) {
1846 "fields_name.size() != field_ents.size()");
1848 for (
auto &f : fields_name) {
1850 if (*it_field_ents) {
1851 for (
auto p = (*it_field_ents)->pair_begin();
1852 p != (*it_field_ents)->pair_end(); ++p) {
1853 a00_fields_uids.emplace_back(
1864 std::vector<const DiagBlockIndex::Indexes *> fe_blocks;
1865 std::vector<const DiagBlockIndex::Indexes *> a00_blocks;
1866 std::vector<const DiagBlockIndex::Indexes *> a01_blocks;
1867 std::vector<const DiagBlockIndex::Indexes *> a10_blocks;
1868 std::vector<const DiagBlockIndex::Indexes *> a11_blocks;
1874 std::vector<int> glob_idx;
1875 std::vector<int> row_glob_idx, col_glob_idx;
1877 auto &block_index = ctx->
blockIndex.get<1>();
1878 for (
auto it = block_index.begin(); it != block_index.end();) {
1887 auto last_uid = it->getFEUId();
1888 while (it != block_index.end() && it->getFEUId() == last_uid) {
1889 fe_blocks.push_back(&*it);
1893 for (
auto b_ptr : fe_blocks) {
1894 auto check_id = [&](
auto uid) {
1895 for (
auto &uid_pair : a00_fields_uids) {
1896 if (uid >= uid_pair.first && uid < uid_pair.second) {
1902 auto r00 = check_id(b_ptr->getRowUId());
1903 auto c00 = check_id(b_ptr->getColUId());
1905 a00_blocks.push_back(b_ptr);
1906 }
else if (r00 && !c00) {
1907 a01_blocks.push_back(b_ptr);
1908 }
else if (!r00 && c00) {
1909 a10_blocks.push_back(b_ptr);
1911 a11_blocks.push_back(b_ptr);
1915 for (
auto b : a11_blocks) {
1917 row_glob_idx.resize(b->getNbRows());
1918 col_glob_idx.resize(b->getNbCols());
1920 std::iota(row_glob_idx.begin(), row_glob_idx.end(), b->getRow());
1921 std::iota(col_glob_idx.begin(), col_glob_idx.end(), b->getCol());
1922 CHKERR AOApplicationToPetsc(ao, col_glob_idx.size(), col_glob_idx.data());
1923 CHKERR AOApplicationToPetsc(ao, row_glob_idx.size(), row_glob_idx.data());
1925 auto shift = b->getMatShift();
1929 col_glob_idx.size(), col_glob_idx.data(), ptr,
1934 col_glob_idx.size(), col_glob_idx.data(), ptr,
1942 <<
"a00_blocks.size() " << a00_blocks.size();
1944 <<
"a01_blocks.size() " << a01_blocks.size();
1946 <<
"a10_blocks.size() " << a10_blocks.size();
1948 <<
"a11_blocks.size() " << a11_blocks.size();
1951 if(a00_blocks.size()) {
1953 for (
auto r : a00_blocks) {
1954 auto r_uid = r->getRowUId();
1955 auto range = ctx->
blockIndex.get<2>().equal_range(r_uid);
1956 for (
auto it = range.first; it != range.second; ++it) {
1957 if (it->getFEUId() != last_uid && it->getColUId() != r_uid) {
1958 a01_blocks.push_back(&*it);
1963 for (
auto c : a00_blocks) {
1964 auto c_uid =
c->getColUId();
1965 auto range = ctx->
blockIndex.get<3>().equal_range(c_uid);
1966 for (
auto it = range.first; it != range.second; ++it) {
1967 if (it->getFEUId() != last_uid && it->getRowUId() != c_uid) {
1968 a10_blocks.push_back(&*it);
1973 auto sort = [](
auto &blocks) {
1974 std::sort(blocks.begin(), blocks.end(), [](
auto a,
auto b) {
1975 if (a->getRowUId() == b->getRowUId())
1976 return a->getColUId() < b->getColUId();
1978 return a->getRowUId() < b->getRowUId();
1988 <<
"a01_blocks.size() " << a01_blocks.size();
1990 <<
"a10_blocks.size() " << a10_blocks.size();
1994 auto set_local_indices = [](
auto &row_blocks,
auto &col_blocks) {
1996 std::map<UId, std::tuple<int, int, int>>
1998 for (
auto b : row_blocks) {
1999 if (block_indexing.find(b->getRowUId()) == block_indexing.end()) {
2000 block_indexing[b->getRowUId()] =
2001 std::make_tuple(b->getNbRows(), b->getNbRows(), b->getRow());
2004 for (
auto b : col_blocks) {
2005 if (block_indexing.find(b->getColUId()) == block_indexing.end()) {
2006 block_indexing[b->getColUId()] =
2007 std::make_tuple(b->getNbCols(), b->getNbCols(), b->getCol());
2011 int mat_block_size = 0;
2012 for (
auto &b : block_indexing) {
2014 auto &[index, size, glob] = b.second;
2015 index = mat_block_size;
2016 mat_block_size += size;
2018 return std::make_pair(block_indexing, mat_block_size);
2021 auto set_block = [&](
auto &blocks,
auto &row_indexing,
auto &col_indexing,
2023 int col_block_size) {
2024 std::vector<std::tuple<int, int, int, int, int>> block_data;
2025 block_data.reserve(blocks.size());
2026 for (
auto s : blocks) {
2027 auto ruid = s->getRowUId();
2028 auto cuid = s->getColUId();
2029 auto &rbi = row_indexing.at(ruid);
2030 auto &cbi = col_indexing.at(cuid);
2031 if (
auto shift = s->getMatShift(); shift != -1) {
2032 block_data.push_back(std::make_tuple(
2036 s->getNbRows(), s->getNbCols(),
2038 get<0>(rbi), get<0>(cbi))
2043 block_mat.resize(row_block_size, col_block_size,
false);
2045 for (
auto &bd : block_data) {
2046 auto &[shift, nb_rows, nb_cols, ridx, cidx] = bd;
2048 for (
auto i = 0;
i != nb_rows; ++
i, ptr += nb_cols) {
2050 cblas_dcopy(nb_cols, ptr, 1, sub_ptr, 1);
2054 for (
auto &bd : block_data) {
2055 auto &[shift, nb_rows, nb_cols, ridx, cidx] = bd;
2057 for (
auto i = 0;
i != nb_rows; ++
i, ptr += nb_cols) {
2059 cblas_daxpy(nb_cols, 1., ptr, 1, sub_ptr, 1);
2065 auto [a00_indexing, a00_size] = set_local_indices(a00_blocks, a00_blocks);
2066 auto [a11_indexing, a11_size] = set_local_indices(a10_blocks, a01_blocks);
2070 <<
"a00_indexing.size() " << a00_indexing.size() <<
" a00_size "
2073 <<
"a11_indexing.size() " << a11_indexing.size() <<
" a11_size "
2077 set_block(a00_blocks, a00_indexing, a00_indexing, block_mat_a00, a00_size,
2079 set_block(a01_blocks, a00_indexing, a11_indexing, block_mat_a01, a00_size,
2081 set_block(a10_blocks, a11_indexing, a00_indexing, block_mat_a10, a11_size,
2083 block_mat_a11.resize(a11_size, a11_size,
false);
2084 block_mat_a11.clear();
2086 block_mat_a00 = trans(block_mat_a00);
2087 block_mat_a01 = trans(block_mat_a01);
2089 ipiv.resize(block_mat_a00.size1(),
false);
2091 lapack_dgesv(block_mat_a00.size1(), block_mat_a01.size1(),
2092 block_mat_a00.data().data(), block_mat_a00.size2(),
2093 &*ipiv.data().begin(), block_mat_a01.data().data(),
2094 block_mat_a01.size2());
2097 "lapack error info = %d", info);
2100 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
2101 block_mat_a11.size1(), block_mat_a11.size2(),
2102 block_mat_a00.size1(), -1.,
2104 block_mat_a10.data().data(), block_mat_a10.size2(),
2106 block_mat_a01.data().data(), block_mat_a01.size2(),
2108 0., block_mat_a11.data().data(), block_mat_a11.size2()
2113 std::vector<int> glob_idx(block_mat_a11.size1());
2114 for (
auto &r : a11_indexing) {
2115 auto [r_index, r_size, r_glob] = r.second;
2116 std::iota(&glob_idx[idx], &glob_idx[idx + r_size], r_glob);
2119 CHKERR AOApplicationToPetsc(ao, glob_idx.size(), glob_idx.data());
2121 glob_idx.data(), block_mat_a11.data().data(),
2134 CHKERR MatShellGetContext(
B, (
void **)&ctx);
2138boost::shared_ptr<std::vector<double>>
2141 CHKERR MatShellGetContext(
B, (
void **)&ctx);
2150 boost::shared_ptr<std::vector<double>> data_blocks_ptr) {
2174 auto get_row_data = [&]() -> std::pair<bool, VectorInt> {
2176 if (
auto stored_data_ptr =
2178 MOFEM_LOG(
"SELF", Sev::warning) <<
"Can lead to unhandled behaviour";
2179 return std::make_pair(
true, stored_data_ptr->entityIndices);
2182 return std::make_pair(
false, row_data.
getIndices());
2185 auto row_indices = get_row_data();
2187 std::vector<int> ent_row_indices;
2188 std::vector<int> ent_col_indices;
2191 if (
auto r_cache = rent->entityCacheRowDofs.lock()) {
2193 auto &row_uid = rent->getLocalUniqueId();
2194 auto &row_ind = row_indices.second;
2197 if (
auto c_cache = cent->entityCacheColDofs.lock()) {
2199 auto &col_uid = cent->getLocalUniqueId();
2204 boost::make_tuple(row_uid, col_uid)
2213 <<
"missing block: row "
2217 "Block not allocated");
2222 auto ri = row_ind.begin();
2223 auto rie = row_ind.end();
2225 auto shift = shift_extractor(&*it);
2226 auto s_mat = &(*data_blocks_ptr)[shift];
2228 auto get_ent_indices = [](
auto &cache,
auto &ind) {
2230 ind.reserve(std::distance(cache->loHi[0], cache->loHi[1]));
2231 for (
auto e = cache->loHi[0]; e != cache->loHi[1]; ++e) {
2232 auto glob = (*e)->getPetscGlobalDofIdx();
2234 ind.push_back(glob);
2238 get_ent_indices(r_cache, ent_row_indices);
2239 if (ent_row_indices.empty())
2241 get_ent_indices(c_cache, ent_col_indices);
2242 if (ent_col_indices.empty())
2245 if (mat.size1() == ent_row_indices.size() &&
2246 mat.size2() == ent_col_indices.size()) {
2248 if (iora == ADD_VALUES) {
2249 cblas_daxpy(mat.data().size(), 1.0, mat.data().data(), 1, s_mat,
2252 cblas_dcopy(mat.data().size(), mat.data().data(), 1, s_mat, 1);
2258 for (
auto re : ent_row_indices) {
2259 ri = std::find(ri, rie, re);
2260 if (!(ri == rie && *ri != -1)) {
2262 auto ci = col_ind.begin();
2263 auto cie = col_ind.end();
2266 for (
auto ce : ent_col_indices) {
2267 ci = std::find(ci, cie, ce);
2268 if (!(ci == cie && *ci != -1)) {
2269 auto &
m = s_mat[row * ent_col_indices.size() + col];
2270 if (iora == ADD_VALUES) {
2271 m += mat(std::distance(row_ind.begin(), ri),
2272 std::distance(col_ind.begin(), ci));
2274 m = mat(std::distance(row_ind.begin(), ri),
2275 std::distance(col_ind.begin(), ci));
2303 PetscBool is_mat_shell = PETSC_FALSE;
2304 PetscObjectTypeCompare((PetscObject)M, MATSHELL, &is_mat_shell);
2307 CHKERR MatShellGetContext(M, (
void **)&ctx);
2309 ctx, row_data, col_data, mat, iora,
2324 PetscBool is_mat_shell = PETSC_FALSE;
2325 PetscObjectTypeCompare((PetscObject)M, MATSHELL, &is_mat_shell);
2328 CHKERR MatShellGetContext(M, (
void **)&ctx);
2331 "No preconditionerBlocksPtr");
2333 ctx, row_data, col_data, mat, iora,
2346 boost::shared_ptr<BlockStructure> block_mat_data_ptr,
2348 std::vector<std::string> fields_names,
2349 std::vector<boost::shared_ptr<Range>>
2351 bool add_preconditioner_block) {
2353 if (!block_mat_data_ptr) {
2357 if (fields_names.size() != field_ents.size())
2359 "fields_names.size() != field_ents.size()");
2361 auto [schur_dm, block_dm] = dms;
2366 auto schur_dofs_row = schur_prb->getNumeredRowDofsPtr();
2367 auto schur_dofs_col = schur_prb->getNumeredColDofsPtr();
2368 auto block_dofs_row = block_prb->getNumeredRowDofsPtr();
2369 auto block_dofs_col = block_prb->getNumeredColDofsPtr();
2371 auto ao_schur_row = schur_prb->getSubData()->getSmartRowMap();
2372 auto ao_schur_col = schur_prb->getSubData()->getSmartColMap();
2373 auto ao_block_row = block_prb->getSubData()->getSmartRowMap();
2374 auto ao_block_col = block_prb->getSubData()->getSmartColMap();
2380 CHKERR VecSetDM(schur_vec_x, PETSC_NULLPTR);
2381 CHKERR VecSetDM(block_vec_x, PETSC_NULLPTR);
2382 CHKERR VecSetDM(schur_vec_y, PETSC_NULLPTR);
2383 CHKERR VecSetDM(block_vec_y, PETSC_NULLPTR);
2385 auto find_field_ent = [&](
auto uid,
auto prb,
auto rc) {
2386 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
2390 dofs = prb->getNumeredRowDofsPtr();
2393 dofs = prb->getNumeredColDofsPtr();
2402 return boost::shared_ptr<NumeredDofEntity>();
2408 return boost::shared_ptr<NumeredDofEntity>();
2411 std::array<boost::shared_ptr<BlockStructure>, 4> data_ptrs;
2413 for (
auto r = 0; r != 3; ++r) {
2414 data_ptrs[r] = boost::make_shared<BlockStructure>();
2415 data_ptrs[r]->dataBlocksPtr = block_mat_data_ptr->dataBlocksPtr;
2417 data_ptrs[3] = boost::make_shared<BlockStructure>();
2418 data_ptrs[3]->dataBlocksPtr = block_mat_data_ptr->dataBlocksPtr;
2420 data_ptrs[0]->ghostX = schur_vec_x;
2421 data_ptrs[0]->ghostY = schur_vec_y;
2422 data_ptrs[1]->ghostX = block_vec_x;
2423 data_ptrs[1]->ghostY = schur_vec_y;
2424 data_ptrs[2]->ghostX = schur_vec_x;
2425 data_ptrs[2]->ghostY = block_vec_y;
2426 data_ptrs[3]->ghostX = block_vec_x;
2427 data_ptrs[3]->ghostY = block_vec_y;
2430 for (
auto &d : block_mat_data_ptr->blockIndex.get<0>()) {
2432 auto insert = [&](
auto &
m,
auto &dof_r,
auto &dof_c,
auto &d) {
2436 d.getRowUId(), d.getColUId(), d.getFEUId(),
2438 dof_r->getPetscGlobalDofIdx(), dof_c->getPetscGlobalDofIdx(),
2440 d.getNbRows(), d.getNbCols(),
2442 dof_r->getPetscLocalDofIdx(), dof_c->getPetscLocalDofIdx(),
2449 auto dof_schur_row = find_field_ent(d.getRowUId(), schur_prb,
ROW);
2450 auto dof_schur_col = find_field_ent(d.getColUId(), schur_prb,
COL);
2451 auto dof_block_row = find_field_ent(d.getRowUId(), block_prb,
ROW);
2452 auto dof_block_col = find_field_ent(d.getColUId(), block_prb,
COL);
2454 if (dof_schur_row && dof_schur_col) {
2455 insert(data_ptrs[0]->blockIndex, dof_schur_row, dof_schur_col, d);
2458 if (dof_schur_row && dof_block_col) {
2459 insert(data_ptrs[1]->blockIndex, dof_schur_row, dof_block_col, d);
2462 if (dof_block_row && dof_schur_col) {
2463 insert(data_ptrs[2]->blockIndex, dof_block_row, dof_schur_col, d);
2466 if (dof_block_row && dof_block_col) {
2467 insert(data_ptrs[3]->blockIndex, dof_block_row, dof_block_col, d);
2474 auto set_up_a00_data = [&](
auto inv_block_data) {
2477 if (add_preconditioner_block) {
2478 auto preconditioned_block = boost::make_shared<std::vector<double>>(
2479 inv_block_data->dataBlocksPtr->size(), 0);
2480 inv_block_data->preconditionerBlocksPtr = preconditioned_block;
2481 inv_block_data->multiplyByPreconditioner =
true;
2482 block_mat_data_ptr->preconditionerBlocksPtr =
2483 inv_block_data->preconditionerBlocksPtr;
2484 block_mat_data_ptr->multiplyByPreconditioner =
false;
2490 CHKERR set_up_a00_data(data_ptrs[3]);
2493 CHKERR PetscObjectGetComm((PetscObject)schur_dm, &comm);
2495 auto create_shell_mat = [&](
auto nb_r_loc,
auto nb_c_loc,
auto nb_r_glob,
2496 auto nb_c_glob,
auto data_ptr) {
2498 CHKERR MatCreateShell(comm, nb_r_loc, nb_c_loc, nb_r_glob, nb_c_glob,
2499 (
void *)data_ptr.get(), &mat_raw);
2504 auto schur_nb_global = schur_prb->getNbDofsRow();
2505 auto block_nb_global = block_prb->getNbDofsRow();
2506 auto schur_nb_local = schur_prb->getNbLocalDofsRow();
2507 auto block_nb_local = block_prb->getNbLocalDofsRow();
2509 std::array<SmartPetscObj<Mat>, 4> mats_array;
2511 create_shell_mat(schur_nb_local, schur_nb_local, schur_nb_global,
2512 schur_nb_global, data_ptrs[0]);
2514 create_shell_mat(schur_nb_local, block_nb_local, schur_nb_global,
2515 block_nb_global, data_ptrs[1]);
2517 create_shell_mat(block_nb_local, schur_nb_local, block_nb_global,
2518 schur_nb_global, data_ptrs[2]);
2520 create_shell_mat(block_nb_local, block_nb_local, block_nb_global,
2521 block_nb_global, data_ptrs[3]);
2524 <<
"(0, 0) " << schur_nb_local <<
" " << schur_nb_global <<
" "
2525 << data_ptrs[0]->blockIndex.size();
2527 <<
"(0, 1) " << schur_nb_local <<
" " << block_nb_local <<
" "
2528 << schur_nb_global <<
" " << block_nb_global <<
" "
2529 << data_ptrs[1]->blockIndex.size();
2531 <<
"(1, 0) " << block_nb_local <<
" " << schur_nb_local <<
" "
2532 << block_nb_global <<
" " << schur_nb_global <<
" "
2533 << data_ptrs[2]->blockIndex.size();
2535 <<
"(1, 1) " << block_nb_local <<
" " << block_nb_global <<
" "
2536 << data_ptrs[3]->blockIndex.size();
2540 auto schur_is = schur_prb->getSubData()->getSmartRowIs();
2541 auto block_is = block_prb->getSubData()->getSmartRowIs();
2543 return boost::make_shared<NestSchurData>(
2545 mats_array, data_ptrs, block_mat_data_ptr,
2546 std::make_pair(schur_is, block_is)
2551std::pair<SmartPetscObj<Mat>, boost::shared_ptr<NestSchurData>>
2554 if (!schur_net_data_ptr)
2557 auto [mat_arrays, data_ptrs, block_mat_data_ptr, iss] = *schur_net_data_ptr;
2558 auto [schur_is, block_is] = iss;
2560 std::array<IS, 2> is_a = {schur_is, block_is};
2561 std::array<Mat, 4> mats_a = {
2563 mat_arrays[0], mat_arrays[1], mat_arrays[2], mat_arrays[3]
2568 CHKERR PetscObjectGetComm((PetscObject)mat_arrays[0], &comm);
2573 comm, 2, is_a.data(), 2, is_a.data(), mats_a.data(), &mat_raw
2584OpSchurAssembleBase *
2586 std::vector<boost::shared_ptr<Range>> field_ents,
2588 bool sym_schur,
bool symm_op) {
2591 schur, sym_schur, symm_op);
2594 schur, sym_schur, symm_op);
2598 boost::shared_ptr<std::vector<unsigned char>> marker_ptr,
double diag_val) {
2602OpSchurAssembleBase *
2604 std::vector<boost::shared_ptr<Range>> field_ents,
2607 std::vector<bool> sym_schur,
bool symm_op,
2608 boost::shared_ptr<BlockStructure> diag_blocks) {
2610 fields_name, field_ents, sequence_of_aos.back(), sequence_of_mats.back(),
2611 sym_schur.back(), symm_op);
2614OpSchurAssembleBase *
2616 std::vector<boost::shared_ptr<Range>> field_ents,
2619 std::vector<bool> sym_schur,
2620 std::vector<double> diag_eps,
bool symm_op,
2621 boost::shared_ptr<BlockStructure> diag_blocks) {
2623 fields_name, field_ents, sequence_of_aos.back(), sequence_of_mats.back(),
2624 sym_schur.back(), symm_op);
2630 auto apply = [](PC pc, Vec f, Vec x) {
2633 CHKERR PCGetOperators(pc, &A, &P);
2634 CHKERR MatSolve(P, f, x);
2638 CHKERR PCSetType(pc, PCSHELL);
2639 CHKERR PCShellSetApply(pc, apply);
2640 CHKERR PCShellSetName(pc,
"MoFEMSchurBlockPC");
2723 col_data, mat, iora);
2755 M, row_data, col_data, mat, iora);
2771 if (block_mat_data->multiplyByPreconditioner) {
2772 block_mat_data->multiplyByPreconditioner =
false;
2774 if (!block_mat_data->preconditionerBlocksPtr) {
2776 "preconditionerBlocksPtr not set");
2778 block_mat_data->multiplyByPreconditioner =
true;
2785 std::string filename) {
2789 moab::Interface &moab = core;
2791 ReadUtilIface *iface;
2792 CHKERR moab.query_interface(iface);
2794 auto size = block_mat_data->blockIndex.size();
2797 vector<double *> arrays_coord;
2798 CHKERR iface->get_node_coords(3, 4 * size, 0, starting_vertex, arrays_coord);
2801 CHKERR iface->get_element_connect(size, 4, MBQUAD, 0, starting_handle, array);
2803 double def_val_nrm2 = 0;
2805 CHKERR moab.tag_get_handle(
"nrm2", 1, MB_TYPE_DOUBLE, th_nrm2,
2806 MB_TAG_CREAT | MB_TAG_DENSE, &def_val_nrm2);
2808 int def_val_mat_shift = 0;
2810 CHKERR moab.tag_get_handle(
"mat_shift", 1, MB_TYPE_INTEGER, th_mat_shift,
2811 MB_TAG_CREAT | MB_TAG_DENSE, &def_val_mat_shift);
2814 for (
auto &d : block_mat_data->blockIndex) {
2815 auto row = d.getRow();
2816 auto col = d.getCol();
2817 auto nb_rows = d.getNbRows();
2818 auto nb_cols = d.getNbCols();
2819 auto mat_shift = d.getMatShift();
2822 arrays_coord[0][4 *
i + 0] = row;
2823 arrays_coord[1][4 *
i + 0] = col;
2824 arrays_coord[2][4 *
i + 0] = 0;
2827 arrays_coord[0][4 *
i + 1] = row + nb_rows;
2828 arrays_coord[1][4 *
i + 1] = col;
2829 arrays_coord[2][4 *
i + 1] = 0;
2832 arrays_coord[0][4 *
i + 2] = row + nb_rows;
2833 arrays_coord[1][4 *
i + 2] = col + nb_cols;
2834 arrays_coord[2][4 *
i + 2] = 0;
2837 arrays_coord[0][4 *
i + 3] = row;
2838 arrays_coord[1][4 *
i + 3] = col + nb_cols;
2839 arrays_coord[2][4 *
i + 3] = 0;
2842 array[4 *
i + 0] = starting_vertex + 4 *
i + 0;
2843 array[4 *
i + 1] = starting_vertex + 4 *
i + 1;
2844 array[4 *
i + 2] = starting_vertex + 4 *
i + 2;
2845 array[4 *
i + 3] = starting_vertex + 4 *
i + 3;
2847 auto prt = &(*block_mat_data->dataBlocksPtr)[d.getMatShift()];
2848 auto nrm2 = cblas_dnrm2(nb_rows * nb_cols, prt, 1);
2850 CHKERR moab.tag_set_data(th_nrm2, &ele, 1, &nrm2);
2851 CHKERR moab.tag_set_data(th_mat_shift, &ele, 1, &mat_shift);
2856 CHKERR iface->update_adjacencies(starting_handle, size, 4, array);
2857 CHKERR moab.write_file(filename.c_str(),
"VTK",
"");
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define BITFIELDID_SIZE
max number of fields
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
FTensor::Index< 'j', 3 > j
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
UBlasVector< int > VectorInt
implementation of Data Operators for Forces and Sources
auto type_from_handle(const EntityHandle h)
get type from entity handle
MoFEMErrorCode MatSetValues< SchurElemMats >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
MoFEMErrorCode schur_mat_set_values_wrap(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
static auto cmp_uid_lo(const boost::weak_ptr< FieldEntity > &a, const UId &b)
static PetscErrorCode solve_add(Mat mat, Vec x, Vec y)
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
EntityHandle get_id_for_max_type()
static auto cmp_uid_hi(const UId &b, const boost::weak_ptr< FieldEntity > &a)
MoFEMErrorCode MatSetValues< SchurElemMatsBlock >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
EntityHandle get_id_for_min_type()
auto isAllGather(IS is)
IS All gather.
MoFEMErrorCode schurSwitchPreconditioner(boost::shared_ptr< BlockStructure > block_mat_data)
Switch preconditioner.
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
constexpr bool debug_schur
auto id_from_handle(const EntityHandle h)
MoFEMErrorCode MatSetValues< SchurElemMatsPreconditionedBlock >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
std::pair< SmartPetscObj< Mat >, boost::shared_ptr< BlockStructure > > SchurShellMatData
std::pair< SmartPetscObj< Mat >, boost::shared_ptr< NestSchurData > > createSchurNestedMatrix(boost::shared_ptr< NestSchurData > schur_net_data_ptr)
Create a Mat Diag Blocks object.
MoFEMErrorCode MatSetValues< BlockStructure >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
MoFEMErrorCode shell_block_mat_asmb_wrap(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
auto createISGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
Creates a data structure for an index set containing a list of integers.
MoFEMErrorCode setSchurMatSolvePC(SmartPetscObj< PC > pc)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
boost::shared_ptr< std::vector< double > > getBlockMatPrconditionerStorageMat(Mat B)
Get the Block Storage object.
static MoFEMErrorCode mult_schur_block_shell(Mat mat, Vec x, Vec y, InsertMode iora, boost::function< int(DiagBlockIndex::BlockIndex::nth_index< 0 >::type::iterator)> shift_extractor, boost::shared_ptr< std::vector< double > > data_blocks_ptr, bool multiply_by_preconditioner)
static MoFEMErrorCode solve_schur_block_shell(Mat mat, Vec y, Vec x, InsertMode iora)
static PetscErrorCode zero_rows_columns(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
static PetscErrorCode mult_add(Mat mat, Vec x, Vec y)
MoFEMErrorCode shell_block_mat_asmb_wrap_impl(BlockStructure *ctx, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora, boost::function< int(const DiagBlockIndex::Indexes *)> shift_extractor, boost::shared_ptr< std::vector< double > > data_blocks_ptr)
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
OpSchurAssembleBase * createOpSchurZeroRowsAndCols(boost::shared_ptr< std::vector< unsigned char > > marker_ptr, double diag_val)
Create a Op Schur Zero Rows And Cols object.
static PetscErrorCode mat_zero(Mat m)
MoFEMErrorCode schurSaveBlockMesh(boost::shared_ptr< BlockStructure > block_mat_data, std::string filename)
Save block matrix as a mesh.
constexpr int max_gemv_size
auto getInterfacePtr(DM dm)
Get the Interface Ptr object.
static PetscErrorCode mult(Mat mat, Vec x, Vec y)
std::vector< std::pair< std::string, std::vector< SchurFieldPair > > > SchurFEOpsFEandFields
static MoFEMErrorCode setSchurBlockMatOps(Mat mat_raw)
MoFEMErrorCode shell_block_preconditioner_mat_asmb_wrap(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
auto getProblemPtr(DM dm)
get problem pointer from DM
SchurShellMatData createBlockMat(DM dm, boost::shared_ptr< BlockStructure > data)
Create a Schur Mat object.
MoFEMErrorCode assembleBlockMatSchur(MoFEM::Interface &m_field, Mat B, Mat S, std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao)
Assemble Schur matrix.
OpSchurAssembleBase * createOpSchurAssembleBegin()
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
boost::shared_ptr< std::vector< double > > getBlockMatStorageMat(Mat B)
Get the Block Storage object.
SmartPetscObj< Mat > block_mat
FTensor::Index< 'm', 3 > m
bool multiplyByPreconditioner
boost::shared_ptr< std::vector< double > > preconditionerBlocksPtr
boost::shared_ptr< std::vector< double > > parentBlockStructurePtr
boost::shared_ptr< std::vector< double > > dataBlocksPtr
SmartPetscObj< Vec > ghostX
SmartPetscObj< Vec > ghostY
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
Deprecated interface functions.
Indexes(UId uid_row, UId uid_col, UId uid_fe, int row, int col, int nb_rows, int nb_cols, int loc_row, int loc_col, int mat_shift)
int & getMatShift() const
multi_index_container< Indexes, indexed_by< ordered_unique< composite_key< Indexes, const_mem_fun< Indexes, UId, &Indexes::getRowUId >, const_mem_fun< Indexes, UId, &Indexes::getColUId > > >, ordered_non_unique< const_mem_fun< Indexes, UId, &Indexes::getFEUId > >, ordered_non_unique< const_mem_fun< Indexes, UId, &Indexes::getRowUId > >, ordered_non_unique< const_mem_fun< Indexes, UId, &Indexes::getColUId > > > > BlockIndex
virtual ~DiagBlockIndex()=default
BlockIndex blockIndex
blocks indexes storage
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getUniqueIdCalculate(const DofIdx dof, UId ent_uid)
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorFieldEntities & getFieldEntities() const
get field entities
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
[Storage and set boundary conditions]
structure for User Loop Methods on finite elements
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static auto getHandleFromUniqueId(const UId uid)
Get the Handle From Unique Id.
static auto getFieldBitNumberFromUniqueId(const UId uid)
Get the Field Bit Number From Unique Id.
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
OpSchurAssembleBaseImpl()=delete
MoFEMErrorCode assembleSchurMat(Mat S, const UId &uid_row, VectorInt &row_ind, const UId &uid_col, VectorInt &col_ind, MatrixDouble &m, InsertMode iora)
Assemble Schur complement.
OpSchurAssembleBase()=delete
Clear Schur complement internal data.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Assemble Schur complement (Implementation)
SmartPetscObj< AO > schurAO
MoFEMErrorCode doWorkImpl(int side, EntityType type, EntitiesFieldData::EntData &data)
SmartPetscObj< Mat > schurMat
std::vector< boost::shared_ptr< Range > > fieldEnts
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.
std::vector< std::string > fieldsName
Assemble Schur complement.
OpSchurZeroRowsAndCols(boost::shared_ptr< std::vector< unsigned char > > marker_ptr, double diag_val)
boost::shared_ptr< std::vector< unsigned char > > markerPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static MatSetValuesPtr matSetValuesPreconditionedBlockPtr
static MatSetValuesPtr matSetValuesPtr
backend assembly function
boost::function< MoFEMErrorCode( Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)> MatSetValuesPtr
static MatSetValuesPtr matSetValuesBlockPtr
backend assembly block mat function
static auto invertMat(MatrixDouble &m, MatrixDouble &inv)
static auto invertMat(MatrixDouble &m, MatrixDouble &inv)
static MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
static MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Schur complement data storage.
static MoFEMErrorCode assembleStorage(const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
static boost::ptr_vector< VectorInt > colIndices
static boost::ptr_vector< MatrixDouble > locMats
SchurElemMats(const size_t idx, const UId uid_row, const UId uid_col)
virtual ~SchurElemMats()=default
static size_t maxIndexCounter
static boost::ptr_vector< VectorInt > rowIndices
static MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
static boost::ptr_vector< SchurElemMats > schurElemMats
static SchurElemStorage schurL2Storage
multi_index_container< const SchurElemMats *, indexed_by< ordered_unique< tag< uid_mi_tag >, composite_key< SchurElemMats, member< SchurElemMats, const UId, &SchurElemMats::uidRow >, member< SchurElemMats, const UId, &SchurElemMats::uidCol > > >, ordered_non_unique< tag< col_mi_tag >, member< SchurElemMats, const UId, &SchurElemMats::uidCol > > > > SchurElemStorage
static PetscLogEvent MOFEM_EVENT_BlockStructureSetValues
static PetscLogEvent MOFEM_EVENT_BlockStructureMult
static PetscLogEvent MOFEM_EVENT_zeroRowsAndCols
static PetscLogEvent MOFEM_EVENT_opSchurAssembleEnd
static PetscLogEvent MOFEM_EVENT_BlockStructureSolve
static PetscLogEvent MOFEM_EVENT_schurMatSetValues
static PetscLogEvent MOFEM_EVENT_AssembleSchurMat
intrusive_ptr for managing petsc objects