v0.15.0
Loading...
Searching...
No Matches
MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE > Struct Template Reference

Assemble Schur complement (Implementation) More...

Inheritance diagram for MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >:
[legend]
Collaboration diagram for MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >:
[legend]

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

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_namelist of fields
field_entslist of entities on which schur complement is applied (can be empty)
schur_aomap indices to schur matrix indices
schur_matschur matrix
sym_schurtrue if schur is symmetric
symm_optrue if block diagonal is symmetric

Definition at line 583 of file Schur.cpp.

588 : OP(NOSPACE, OP::OPSPACE, symm_op), fieldsName(fields_name),
589 fieldEnts(field_ents), schurAO(schur_ao), schurMat(schur_mat),
590 symSchur(sym_schur) {}
@ NOSPACE
Definition definitions.h:83
SmartPetscObj< AO > schurAO
Definition Schur.cpp:122
OP_SCHUR_ASSEMBLE_BASE OP
Definition Schur.cpp:90
SmartPetscObj< Mat > schurMat
Definition Schur.cpp:123
std::vector< boost::shared_ptr< Range > > fieldEnts
Definition Schur.cpp:121
std::vector< std::string > fieldsName
Definition Schur.cpp:120

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.

595 {
596
598
599#ifndef NDEBUG
600 PetscLogEventBegin(SchurEvents::MOFEM_EVENT_opSchurAssembleEnd, 0, 0, 0, 0);
601#endif
602
603#ifndef NDEBUG
604 if constexpr (debug_schur)
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;
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];
618 if (row_ents) {
619 for (auto p = row_ents->pair_begin(); p != row_ents->pair_end(); ++p) {
620 auto lo_uid =
621 FieldEntity::getLoLocalEntityBitNumber(field_bit, p->first);
622 auto hi_uid =
623 FieldEntity::getHiLocalEntityBitNumber(field_bit, p->second);
624 a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
625 }
626 } else {
628 field_bit, get_id_for_min_type<MBVERTEX>());
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) {
639 return OP::getPtrFE()->mField.get_field_name(field_bit_from_bit_number(
641 };
642 auto list_storage = [&](auto &storage) {
644 int i = 0;
645 for (auto &p : storage) {
646 MOFEM_LOG("SELF", Sev::noisy)
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();
650 ++i;
651 }
653 };
654#endif // NDEBUG
655
656 auto assemble_dense_blocks = [&]() {
657 using matrix_range = ublas::matrix_range<MatrixDouble>;
658 using range = ublas::range;
660 auto &storage = SchurElemMats::schurL2Storage;
661
662 auto assemble_schur = [this](auto &m, auto &uid_row, auto &uid_col,
663 auto *row_ind_ptr, auto *col_ind_ptr) {
665
666 if (schurMat) {
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);
678 MOFEM_LOG_CHANNEL("SELF");
679 if (row_ent_it != field_ents->end())
680 MOFEM_LOG("SELF", Sev::error)
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())
685 MOFEM_LOG("SELF", Sev::error)
686 << "Assemble col entity: " << (*col_ent_it)->getName() << " "
687 << (*col_ent_it)->getEntTypeName() << " side "
688 << (*col_ent_it)->getSideNumber();
689#endif // NDEBUG
690 CHK_THROW_MESSAGE(ierr, "MatSetValues");
691 }
692 }
693
695 };
696
697 auto a00_uids = get_a00_uids();
698
699 auto get_block_indexing = [&](auto &a00_uids) {
700 // iterate over a00 uids and find blocks
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 // create block indexes map for blockMat
725 std::map<UId, std::pair<size_t, const VectorInt *>>
726 block_indexing; // uid block map
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 // set indexes to block
735 int mat_block_size = 0; // size of block matrix
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
766 if (schurAO) {
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());
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);
784 blockMat.clear();
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())
794 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong rlo_uid");
795 if (block_indexing.find(s->uidCol) == block_indexing.end())
796 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong clo_uid");
797#endif // NDEBUG
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));
803 sub_mat = m;
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) {
833 for (auto j = 0; j != blockMat.size2(); ++j) {
834 blockMat(i, j) = 0;
835 }
836 }
837 for (auto j : zero_cols) {
838 for (auto i = 0; i != blockMat.size1(); ++i) {
839 blockMat(i, j) = 0;
840 }
841 }
842 for (auto i : zero_rows) {
843 blockMat(i, i) = 1;
844 }
845
846 CHKERR I::invertMat(blockMat, invMat);
847
848 // clear storage and block list from a00 blocks, no more needed
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
856 if (schurAO) {
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());
861 }
862 if (block_indexing.find(s->uidCol) == block_indexing.end()) {
863 CHKERR AOApplicationToPetsc(schurAO, s->getColInd().size(),
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(
916 boost::make_tuple(cm->first, FieldEntity::getHiBitNumberUId(
917 BITFIELDID_SIZE - 1)));
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())
927 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
928 "Wrong a_lo->uidRow");
929#endif
930
931 auto &a = (*a_lo)->getMat();
932 abM.resize(a.size1(), bM.size2(), false);
933 abM.clear();
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())
948 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
949 "Wrong a_lo->uidRow");
950#endif
951
952 auto &c = (*c_lo)->getMat();
953
954 abcM.resize(abM.size1(), c.size2(), false);
955 abcM.clear();
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) {
971 abcM = trans(abcM);
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 // Assemble Schur complements
987 CHKERR assemble_dense_blocks();
988
989#ifndef NDEBUG
990 if constexpr (debug_schur)
991 MOFEM_LOG("SELF", Sev::noisy) << "Schur assemble done";
992#endif
993
994#ifndef NDEBUG
995 PetscLogEventEnd(SchurEvents::MOFEM_EVENT_opSchurAssembleEnd, 0, 0, 0, 0);
996#endif
997
999}
constexpr double a
#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
Definition definitions.h:31
#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
Definition Schur.cpp:12
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
Definition Schur.cpp:218
static PetscLogEvent MOFEM_EVENT_opSchurAssembleEnd
Definition Schur.hpp:28

Member Data Documentation

◆ abcM

Definition at line 128 of file Schur.cpp.

◆ abM

Definition at line 128 of file Schur.cpp.

◆ blockMat

Definition at line 126 of file Schur.cpp.

◆ bM

Definition at line 128 of file Schur.cpp.

◆ fieldEnts

template<typename OP_SCHUR_ASSEMBLE_BASE >
std::vector<boost::shared_ptr<Range> > MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >::fieldEnts
protected

Definition at line 121 of file Schur.cpp.

◆ fieldsName

template<typename OP_SCHUR_ASSEMBLE_BASE >
std::vector<std::string> MoFEM::OpSchurAssembleEndImpl< OP_SCHUR_ASSEMBLE_BASE >::fieldsName
protected

Definition at line 120 of file Schur.cpp.

◆ invMat

Definition at line 127 of file Schur.cpp.

◆ schurAO

Definition at line 122 of file Schur.cpp.

◆ schurMat

Definition at line 123 of file Schur.cpp.

◆ symSchur

Definition at line 124 of file Schur.cpp.


The documentation for this struct was generated from the following file: