595 {
596
598
599#ifndef NDEBUG
601#endif
602
603#ifndef NDEBUG
605 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble begin -> end";
606#endif
607
608 auto get_a00_uids = [&]() {
609 auto get_field_bit = [&](auto &name) {
610 return OP::getPtrFE()->mField.get_field_bit_number(name);
611 };
612
613 std::vector<std::pair<UId, UId>> a00_uids;
615 for (
auto ss = 0; ss !=
fieldsName.size(); ++ss) {
616 auto field_bit = get_field_bit(
fieldsName[ss]);
618 if (row_ents) {
619 for (auto p = row_ents->pair_begin(); p != row_ents->pair_end(); ++p) {
620 auto lo_uid =
622 auto hi_uid =
624 a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
625 }
626 } else {
631 a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
632 }
633 }
634 return a00_uids;
635 };
636
637#ifndef NDEBUG
638 auto get_field_name = [&](auto uid) {
641 };
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();
651 }
653 };
654#endif
655
656 auto assemble_dense_blocks = [&]() {
657 using matrix_range = ublas::matrix_range<MatrixDouble>;
658 using range = ublas::range;
661
662 auto assemble_schur = [
this](
auto &
m,
auto &uid_row,
auto &uid_col,
663 auto *row_ind_ptr, auto *col_ind_ptr) {
665
667
668 if (
auto ierr = this->assembleSchurMat(
669
670 schurMat, uid_row, *row_ind_ptr, uid_col, *col_ind_ptr,
m,
671 ADD_VALUES
672
673 )) {
674#ifndef NDEBUG
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();
689#endif
691 }
692 }
693
695 };
696
697 auto a00_uids = get_a00_uids();
698
699 auto get_block_indexing = [&](auto &a00_uids) {
700
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;
707
708 auto it =
709 storage.template get<SchurElemMats::uid_mi_tag>().lower_bound(
710 boost::make_tuple(rlo_uid, clo_uid));
711 auto hi_it =
712 storage.template get<SchurElemMats::uid_mi_tag>().upper_bound(
713 boost::make_tuple(rhi_uid, chi_uid));
714
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);
719 }
720 }
721 }
722 }
723
724
725 std::map<UId, std::pair<size_t, const VectorInt *>>
726 block_indexing;
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()));
731 }
732 }
733
734
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();
743 }
744 }
745
746 return std::make_tuple(block_list, block_indexing, mat_block_size);
747 };
748
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);
756 }
757 }
758 return block_list;
759 };
760
761 auto bi = get_block_indexing(a00_uids);
762 auto &[block_list, block_indexing, block_mat_size] = bi;
763
764 if (block_mat_size == 0) {
765
767 for (auto &s : storage) {
769 &*s->getRowInd().begin());
771 &*s->getColInd().begin());
772 }
773 }
774
775 for (auto &s : storage) {
776 auto &
m = s->getMat();
777 CHKERR assemble_schur(
m, s->uidRow, s->uidCol, &(s->getRowInd()),
778 &(s->getColInd()));
779 }
781 }
782
783 blockMat.resize(block_mat_size, block_mat_size,
false);
785
786 auto get_range = [](auto &bi) {
787 return range(bi.first, bi.first + bi.second->size());
788 };
789
790 for (auto &s : block_list) {
791 auto &
m = s->getMat();
792#ifndef NDEBUG
793 if (block_indexing.find(s->uidRow) == block_indexing.end())
795 if (block_indexing.find(s->uidCol) == block_indexing.end())
797#endif
798
799 auto &rbi = block_indexing.at(s->uidRow);
800 auto &cbi = block_indexing.at(s->uidCol);
801
802 auto sub_mat = matrix_range(
blockMat, get_range(rbi), get_range(cbi));
804 }
805
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) {
814 if (*it < 0) {
815 auto idx = bi.first + std::distance(
ind.begin(), it);
816 zeroed_indices.push_back(idx);
817 }
818 }
819 }
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;
824 };
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(); });
831
832 for (
auto i : zero_rows) {
835 }
836 }
837 for (
auto j : zero_cols) {
840 }
841 }
842 for (
auto i : zero_rows) {
844 }
845
847
848
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);
853 }
854 block_list.clear();
855
857 for (auto &s : storage) {
858 if (block_indexing.find(s->uidRow) == block_indexing.end()) {
860 &*s->getRowInd().begin());
861 }
862 if (block_indexing.find(s->uidCol) == block_indexing.end()) {
864 &*s->getColInd().begin());
865 }
866 }
867 }
868
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()),
873 &(s->getColInd()));
874 }
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);
879 }
880 schur_block_list.clear();
881
882 for (
883
884 auto rp_uid_it = a00_uids.begin(); rp_uid_it != a00_uids.end();
885 ++rp_uid_it
886
887 ) {
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;
892
893 auto a_lo_tmp =
894 storage.template get<SchurElemMats::col_mi_tag>().lower_bound(
895 rm->first);
896 auto a_hi =
897 storage.template get<SchurElemMats::col_mi_tag>().upper_bound(
898 rm->first);
899
900 for (
901
902 auto cp_uid_it = (
symSchur) ? rp_uid_it : a00_uids.begin();
903 cp_uid_it != a00_uids.end(); ++cp_uid_it
904
905 ) {
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;
910
911 auto c_lo_tmp =
912 storage.template get<SchurElemMats::uid_mi_tag>().lower_bound(
913 boost::make_tuple(cm->first, 0));
914 auto c_hi =
915 storage.template get<SchurElemMats::uid_mi_tag>().upper_bound(
918
919 auto sub_inv_mat =
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;
923
924 for (auto a_lo = a_lo_tmp; a_lo != a_hi; ++a_lo) {
925#ifndef NDEBUG
926 if (block_indexing.find((*a_lo)->uidRow) != block_indexing.end())
928 "Wrong a_lo->uidRow");
929#endif
930
931 auto &
a = (*a_lo)->getMat();
932 abM.resize(
a.size1(),
bM.size2(),
false);
934 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
935
936 a.size1(),
bM.size2(),
a.size2(), 1.,
937
938 &*
a.data().begin(),
a.size2(),
939
940 &*
bM.data().begin(),
bM.size2(), 0.,
941
942 &*
abM.data().begin(),
abM.size2());
943
944 for (auto c_lo = c_lo_tmp; c_lo != c_hi; ++c_lo) {
945#ifndef NDEBUG
946 if (block_indexing.find((*c_lo)->uidCol) !=
947 block_indexing.end())
949 "Wrong a_lo->uidRow");
950#endif
951
952 auto &
c = (*c_lo)->getMat();
953
954 abcM.resize(
abM.size1(),
c.size2(),
false);
956
957 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
958 abM.size1(),
c.size2(),
abM.size2(), -1.,
959
960 &*
abM.data().begin(),
abM.size2(),
961
962 &*
c.data().begin(),
c.size2(), 0.,
963
964 &*
abcM.data().begin(),
abcM.size2());
965
966 CHKERR assemble_schur(
abcM, (*a_lo)->uidRow, (*c_lo)->uidCol,
967 &((*a_lo)->getRowInd()),
968 &((*c_lo)->getColInd()));
969
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()));
975 }
976 }
977 }
978 }
979 }
980 }
981 }
982
984 };
985
986
987 CHKERR assemble_dense_blocks();
988
989#ifndef NDEBUG
991 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble done";
992#endif
993
994#ifndef NDEBUG
996#endif
997
999}
#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_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#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)
FTensor::Index< 'j', 3 > j
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
EntityHandle get_id_for_max_type()
EntityHandle get_id_for_min_type()
constexpr bool debug_schur
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
FTensor::Index< 'm', 3 > m
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static auto getFieldBitNumberFromUniqueId(const UId uid)
Get the Field Bit Number From Unique Id.
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static SchurElemStorage schurL2Storage
static PetscLogEvent MOFEM_EVENT_opSchurAssembleEnd