6#define MatrixManagerFunctionBegin \
8 MOFEM_LOG_CHANNEL("WORLD"); \
9 MOFEM_LOG_CHANNEL("SYNC"); \
10 MOFEM_LOG_FUNCTION(); \
11 MOFEM_LOG_TAG("SYNC", "MatrixManager"); \
12 MOFEM_LOG_TAG("WORLD", "MatrixManager")
36 MPI_Comm comm = PETSC_COMM_WORLD,
int verbose = 1)
39 typedef FieldEntityEntFiniteElementAdjacencyMap_multiIndex::index<
42 typedef NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type
54 template <
typename TAG>
57 std::vector<PetscInt> &
i, std::vector<PetscInt> &
j,
58 const bool no_diagonals =
true,
int verb =
QUIET)
const;
62 template <
typename TAG>
64 ProblemsByName::iterator p_miit,
66 TAG>::type::iterator mit_row,
67 boost::shared_ptr<FieldEntity> mofem_ent_ptr,
68 std::vector<int> &dofs_col_view,
int verb)
const;
71template <
typename TAG>
73 ProblemsByName::iterator p_miit,
75 TAG>::type::iterator mit_row,
76 boost::shared_ptr<FieldEntity> mofem_ent_ptr,
77 std::vector<int> &dofs_col_view,
int verb)
const {
81 BitRefLevel prb_mask = p_miit->getBitRefLevelMask();
85 dofs_col_view.clear();
87 mofem_ent_ptr->getLocalUniqueId());
88 r.first != r.second; ++r.first) {
90 if (r.first->byWhat &
BYROW) {
92 if ((r.first->entFePtr->getId() & p_miit->getBitFEId()).none()) {
97 const BitRefLevel &fe_bit = r.first->entFePtr->getBitRefLevel();
99 if ((fe_bit & prb_bit).any() && (fe_bit & prb_mask) == fe_bit) {
101 const bool empty_row_block =
102 (empty_field_blocks.first & (*mit_row)->getId()).none();
104 for (
auto &it : r.first->entFePtr->getColFieldEnts()) {
105 if (
auto e = it.lock()) {
106 if (empty_row_block ||
107 (empty_field_blocks.second & e->getId()).none()) {
109 if (
auto cache = e->entityCacheColDofs.lock()) {
110 const auto lo =
cache->loHi[0];
111 const auto hi =
cache->loHi[1];
112 for (
auto vit = lo; vit != hi; ++vit) {
114 const int idx = TAG::get_index(vit);
115 if (PetscLikely(idx >= 0))
116 dofs_col_view.push_back(idx);
119 const DofIdx nb_dofs_col = p_miit->getNbDofsCol();
120 if (PetscUnlikely(idx >= nb_dofs_col)) {
122 <<
"Problem with dof: " << std::endl
123 <<
"Rank " <<
rAnk <<
" : " << *(*vit);
126 "Index of dof larger than number of DOFs in column");
144template <
typename TAG>
146 ProblemsByName::iterator p_miit,
const MatType type,
147 std::vector<PetscInt> &
i, std::vector<PetscInt> &
j,
const bool no_diagonals,
153 auto cache = boost::make_shared<std::vector<EntityCacheNumeredDofs>>(
159 const auto uid = (*it)->getLocalUniqueId();
161 for (
auto lo = r.first; lo != r.second; ++lo) {
163 if ((lo->getBitFEId() & p_miit->getBitFEId()).any()) {
165 const BitRefLevel &prb_bit = p_miit->getBitRefLevel();
166 const BitRefLevel &prb_mask = p_miit->getBitRefLevelMask();
167 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
170 if ((fe_bit & prb_mask) != fe_bit)
172 if ((fe_bit & prb_bit).none())
175 auto dit = p_miit->numeredColDofsPtr->lower_bound(uid);
177 decltype(dit) hi_dit;
178 if (dit != p_miit->numeredColDofsPtr->end())
179 hi_dit = p_miit->numeredColDofsPtr->upper_bound(
184 (*it)->entityCacheColDofs =
185 boost::shared_ptr<EntityCacheNumeredDofs>(
cache, &((*
cache)[idx]));
186 (*cache)[idx].loHi = {dit, hi_dit};
193 using NumeredDofEntitysByIdx =
198 const NumeredDofEntitysByIdx &dofs_row_by_idx =
199 p_miit->numeredRowDofsPtr->get<TAG>();
200 int nb_dofs_row = p_miit->getNbDofsRow();
201 if (nb_dofs_row == 0) {
203 p_miit->getName().c_str());
207 std::map<int, std::vector<int>> adjacent_dofs_on_other_parts;
214 TAG>::type::iterator miit_row,
217 if (TAG::IamNotPartitioned) {
222 CHKERR PetscLayoutSetBlockSize(layout, 1);
223 CHKERR PetscLayoutSetSize(layout, nb_dofs_row);
224 CHKERR PetscLayoutSetUp(layout);
225 PetscInt rstart, rend;
226 CHKERR PetscLayoutGetRange(layout, &rstart, &rend);
227 CHKERR PetscLayoutDestroy(&layout);
231 <<
"row lower " << rstart <<
" row upper " << rend;
235 miit_row = dofs_row_by_idx.lower_bound(rstart);
236 hi_miit_row = dofs_row_by_idx.lower_bound(rend);
237 if (std::distance(miit_row, hi_miit_row) != rend - rstart) {
240 "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
241 "rstart (%d != %d - %d = %d) ",
242 std::distance(miit_row, hi_miit_row), rend, rstart, rend - rstart);
252 std::vector<std::vector<int>> dofs_vec(
sIze);
254 boost::shared_ptr<FieldEntity> mofem_ent_ptr;
255 std::vector<int> dofs_col_view;
258 TAG>::type::iterator mit_row,
261 mit_row = dofs_row_by_idx.begin();
262 hi_mit_row = dofs_row_by_idx.end();
263 for (; mit_row != hi_mit_row; mit_row++) {
271 unsigned char pstatus = (*mit_row)->getPStatus();
272 if ((pstatus & PSTATUS_NOT_OWNED) &&
273 (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED))) {
275 bool get_adj_col =
true;
277 if (mofem_ent_ptr->getLocalUniqueId() ==
278 (*mit_row)->getFieldEntityPtr()->getLocalUniqueId()) {
285 mofem_ent_ptr = (*mit_row)->getFieldEntityPtr();
286 CHKERR getEntityAdjacenies<TAG>(p_miit, mit_row, mofem_ent_ptr,
287 dofs_col_view, verb);
290 std::sort(dofs_col_view.begin(), dofs_col_view.end());
291 std::vector<int>::iterator new_end =
292 std::unique(dofs_col_view.begin(), dofs_col_view.end());
293 int new_size = std::distance(dofs_col_view.begin(), new_end);
294 dofs_col_view.resize(new_size);
298 int owner = (*mit_row)->getOwnerProc();
299 dofs_vec[owner].emplace_back(TAG::get_index(mit_row));
300 dofs_vec[owner].emplace_back(
301 dofs_col_view.size());
303 dofs_vec[owner].insert(dofs_vec[owner].end(), dofs_col_view.begin(),
304 dofs_col_view.end());
310 std::vector<int> dofs_vec_length(
sIze);
311 for (
int proc = 0; proc <
sIze; proc++) {
313 if (!dofs_vec[proc].empty()) {
315 dofs_vec_length[proc] = dofs_vec[proc].size();
320 dofs_vec_length[proc] = 0;
324 std::vector<MPI_Status> status(
sIze);
328 CHKERR PetscGatherNumberOfMessages(comm, NULL, &dofs_vec_length[0],
335 CHKERR PetscGatherMessageLengths(comm, nsends, nrecvs, &dofs_vec_length[0],
340 CHKERR PetscCommGetNewTag(comm, &tag);
345 MPI_Request *r_waits;
350 CHKERR PetscPostIrecvInt(comm, tag, nrecvs, onodes, olengths, &rbuf,
353 MPI_Request *s_waits;
354 CHKERR PetscMalloc1(nsends, &s_waits);
357 for (
int proc = 0, kk = 0; proc <
sIze; proc++) {
358 if (!dofs_vec_length[proc])
360 CHKERR MPI_Isend(&(dofs_vec[proc])[0],
361 dofs_vec_length[proc],
363 tag, comm, s_waits + kk);
369 CHKERR MPI_Waitall(nrecvs, r_waits, &status[0]);
373 CHKERR MPI_Waitall(nsends, s_waits, &status[0]);
376 for (
int kk = 0; kk < nrecvs; kk++) {
378 int len = olengths[kk];
379 int *data_from_proc = rbuf[kk];
381 for (
int ii = 0; ii < len;) {
383 int row_idx = data_from_proc[ii++];
384 int nb_adj_dofs = data_from_proc[ii++];
388 DofByGlobalPetscIndex::iterator dit;
394 "dof %d can not be found in problem", row_idx);
398 for (
int jj = 0; jj < nb_adj_dofs; jj++) {
399 adjacent_dofs_on_other_parts[row_idx].push_back(data_from_proc[ii++]);
405 CHKERR PetscFree(s_waits);
406 CHKERR PetscFree(rbuf[0]);
408 CHKERR PetscFree(r_waits);
410 CHKERR PetscFree(olengths);
412 miit_row = dofs_row_by_idx.begin();
413 hi_miit_row = dofs_row_by_idx.end();
415 CHKERR PetscCommDestroy(&comm);
418 boost::shared_ptr<FieldEntity> mofem_ent_ptr;
419 int row_last_evaluated_idx = -1;
421 std::vector<int> dofs_vec;
422 std::vector<int> dofs_col_view;
425 int nb_loc_row_from_iterators = 0;
426 unsigned int rows_to_fill = p_miit->getNbLocalDofsRow();
427 i.reserve(rows_to_fill + 1);
428 for (; miit_row != hi_miit_row; miit_row++) {
430 if (!TAG::IamNotPartitioned) {
431 if (
static_cast<int>((*miit_row)->getPart()) !=
rAnk)
435 nb_loc_row_from_iterators++;
438 i.push_back(
j.size());
445 : (mofem_ent_ptr->getLocalUniqueId() !=
446 (*miit_row)->getFieldEntityPtr()->getLocalUniqueId())) {
449 mofem_ent_ptr = (*miit_row)->getFieldEntityPtr();
450 CHKERR getEntityAdjacenies<TAG>(p_miit, miit_row, mofem_ent_ptr,
451 dofs_col_view, verb);
452 row_last_evaluated_idx = TAG::get_index(miit_row);
456 dofs_vec.insert(dofs_vec.end(), dofs_col_view.begin(),
457 dofs_col_view.end());
459 unsigned char pstatus = (*miit_row)->getPStatus();
461 std::map<int, std::vector<int>>::iterator mit;
462 mit = adjacent_dofs_on_other_parts.find(row_last_evaluated_idx);
463 if (mit == adjacent_dofs_on_other_parts.end()) {
471 dofs_vec.insert(dofs_vec.end(), mit->second.begin(),
477 sort(dofs_vec.begin(), dofs_vec.end());
478 std::vector<int>::iterator new_end =
479 unique(dofs_vec.begin(), dofs_vec.end());
480 int new_size = std::distance(dofs_vec.begin(), new_end);
481 dofs_vec.resize(new_size);
485 if (
j.capacity() <
j.size() + dofs_vec.size()) {
488 unsigned int nb_nonzero =
j.size() + dofs_vec.size();
489 unsigned int average_row_fill = nb_nonzero /
i.size() + 1;
490 j.reserve(rows_to_fill * average_row_fill);
493 auto hi_diit = dofs_vec.end();
494 for (
auto diit = dofs_vec.begin(); diit != hi_diit; diit++) {
497 if (*diit == TAG::get_index(miit_row)) {
506 i.push_back(
j.size());
508 if (strcmp(type, MATMPIADJ) == 0) {
511 if (
i.size() - 1 != (
unsigned int)nb_loc_row_from_iterators) {
514 "Number of rows from iterator is different than size of rows in "
516 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
518 (
unsigned int)nb_loc_row_from_iterators,
i.size() - 1);
521 }
else if (strcmp(type, MATMPIAIJ) == 0) {
524 if (
i.size() - 1 != (
unsigned int)nb_loc_row_from_iterators) {
527 "Number of rows from iterator is different than size of rows in "
529 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
531 (
unsigned int)nb_loc_row_from_iterators,
i.size() - 1);
533 PetscInt nb_local_dofs_row = p_miit->getNbLocalDofsRow();
534 if ((
unsigned int)nb_local_dofs_row !=
i.size() - 1) {
537 "Number of rows is different than size of rows in compressed row "
538 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
540 (
unsigned int)nb_local_dofs_row,
i.size() - 1);
543 }
else if (strcmp(type, MATAIJ) == 0) {
546 if (
i.size() - 1 != (
unsigned int)nb_loc_row_from_iterators) {
549 "Number of rows form iterator is different than size of rows in "
551 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
553 (
unsigned int)nb_loc_row_from_iterators,
i.size() - 1);
555 PetscInt nb_local_dofs_row = p_miit->getNbLocalDofsRow();
556 if ((
unsigned int)nb_local_dofs_row !=
i.size() - 1) {
559 "Number of rows is different than size of rows in compressed row "
560 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
562 (
unsigned int)nb_local_dofs_row,
i.size() - 1);
565 }
else if (strcmp(type, MATAIJCUSPARSE) == 0) {
567 if (
i.size() - 1 != (
unsigned int)nb_loc_row_from_iterators) {
568 SETERRQ(
get_comm(), PETSC_ERR_ARG_SIZ,
"data inconsistency");
570 PetscInt nb_local_dofs_row = p_miit->getNbLocalDofsRow();
571 if ((
unsigned int)nb_local_dofs_row !=
i.size() - 1) {
572 SETERRQ(
get_comm(), PETSC_ERR_ARG_SIZ,
"data inconsistency");
577 SETERRQ(
get_comm(), PETSC_ERR_ARG_NULL,
"not implemented");
592 : cOre(const_cast<
MoFEM::
Core &>(core)) {
593 PetscLogEventRegister(
"MatrixManagerCreateMPIAIJ", 0,
595 PetscLogEventRegister(
"MatrixManagerCreateMPIAIJWithArrays", 0,
597 PetscLogEventRegister(
"MatrixManagerCreateMPIAdjWithArrays", 0,
599 PetscLogEventRegister(
"MatrixManagerCreateMPIAIJCUSPARSEWithArrays", 0,
601 PetscLogEventRegister(
"MatrixManagerCreateSeqAIJCUSPARSEWithArrays", 0,
603 PetscLogEventRegister(
"MatrixManagerCreateSeqAIJWithArrays", 0,
605 PetscLogEventRegister(
"MatrixManagerCheckMPIAIJWithArraysMatrixFillIn", 0,
611 const std::string name, Mat *Aij,
int verb) {
621 auto p_miit = prb.find(name);
622 if (p_miit == prb.end()) {
624 "problem < %s > is not found (top tip: check spelling)",
628 std::vector<int> i_vec, j_vec;
629 j_vec.reserve(10000);
631 p_miit, MATMPIAIJ, i_vec, j_vec,
false, verb);
633 int nb_row_dofs = p_miit->getNbDofsRow();
634 int nb_col_dofs = p_miit->getNbDofsCol();
635 int nb_local_dofs_row = p_miit->getNbLocalDofsRow();
636 int nb_local_dofs_col = p_miit->getNbLocalDofsCol();
638 CHKERR ::MatCreateMPIAIJWithArrays(
639 m_field.
get_comm(), nb_local_dofs_row, nb_local_dofs_col, nb_row_dofs,
640 nb_col_dofs, &*i_vec.begin(), &*j_vec.begin(), PETSC_NULL, Aij);
648MatrixManager::createMPIAIJCUSPARSEWithArrays<PetscGlobalIdx_mi_tag>(
649 const std::string name, Mat *Aij,
int verb) {
659 auto p_miit = prb.find(name);
660 if (p_miit == prb.end()) {
662 "problem < %s > is not found (top tip: check spelling)",
666 std::vector<int> i_vec, j_vec;
667 j_vec.reserve(10000);
669 p_miit, MATAIJCUSPARSE, i_vec, j_vec,
false, verb);
671 int nb_local_dofs_row = p_miit->getNbLocalDofsRow();
672 int nb_local_dofs_col = p_miit->getNbLocalDofsCol();
674 auto get_layout = [&]() {
675 int start_ranges, end_ranges;
678 CHKERR PetscLayoutSetBlockSize(layout, 1);
679 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs_col);
680 CHKERR PetscLayoutSetUp(layout);
681 CHKERR PetscLayoutGetRange(layout, &start_ranges, &end_ranges);
682 CHKERR PetscLayoutDestroy(&layout);
683 return std::make_pair(start_ranges, end_ranges);
686 auto get_nnz = [&](
auto &d_nnz,
auto &o_nnz) {
688 auto layout = get_layout();
690 for (
int i = 0;
i != nb_local_dofs_row; ++
i) {
691 for (;
j != i_vec[
i + 1]; ++
j) {
692 if (j_vec[
j] < layout.second && j_vec[
j] >= layout.first)
701 std::vector<int> d_nnz(nb_local_dofs_row, 0), o_nnz(nb_local_dofs_row, 0);
702 CHKERR get_nnz(d_nnz, o_nnz);
704#ifdef PETSC_HAVE_CUDA
705 CHKERR ::MatCreateAIJCUSPARSE(m_field.
get_comm(), nb_local_dofs_row,
706 nb_local_dofs_col, nb_row_dofs, nb_col_dofs, 0,
707 &*d_nnz.begin(), 0, &*o_nnz.begin(), Aij);
710 "Error: To use this matrix type compile PETSc with CUDA.");
720MatrixManager::createSeqAIJCUSPARSEWithArrays<PetscLocalIdx_mi_tag>(
721 const std::string name, Mat *Aij,
int verb) {
724#ifdef PETSC_HAVE_CUDA
728 "Not implemented type of matrix yet, try MPI version (aijcusparse)");
736MatrixManager::createMPIAIJ<PetscGlobalIdx_mi_tag>(
const std::string name,
737 Mat *Aij,
int verb) {
746 auto p_miit = prb.find(name);
747 if (p_miit == prb.end()) {
749 "problem < %s > is not found (top tip: check spelling)",
753 std::vector<int> i_vec, j_vec;
754 j_vec.reserve(10000);
756 p_miit, MATMPIAIJ, i_vec, j_vec,
false, verb);
758 int nb_row_dofs = p_miit->getNbDofsRow();
759 int nb_col_dofs = p_miit->getNbDofsCol();
760 int nb_local_dofs_row = p_miit->getNbLocalDofsRow();
761 int nb_local_dofs_col = p_miit->getNbLocalDofsCol();
763 auto get_layout = [&]() {
764 int start_ranges, end_ranges;
767 CHKERR PetscLayoutSetBlockSize(layout, 1);
768 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs_col);
769 CHKERR PetscLayoutSetUp(layout);
770 CHKERR PetscLayoutGetRange(layout, &start_ranges, &end_ranges);
771 CHKERR PetscLayoutDestroy(&layout);
772 return std::make_pair(start_ranges, end_ranges);
775 auto get_nnz = [&](
auto &d_nnz,
auto &o_nnz) {
777 auto layout = get_layout();
779 for (
int i = 0;
i != nb_local_dofs_row; ++
i) {
780 for (;
j != i_vec[
i + 1]; ++
j) {
781 if (j_vec[
j] < layout.second && j_vec[
j] >= layout.first)
790 std::vector<int> d_nnz(nb_local_dofs_row, 0), o_nnz(nb_local_dofs_row, 0);
791 CHKERR get_nnz(d_nnz, o_nnz);
794 CHKERR MatSetSizes(*Aij, nb_local_dofs_row, nb_local_dofs_col, nb_row_dofs,
796 CHKERR MatSetType(*Aij, MATMPIAIJ);
797 CHKERR MatMPIAIJSetPreallocation(*Aij, 0, &*d_nnz.begin(), 0,
806MatrixManager::createMPIAdjWithArrays<Idx_mi_tag>(
const std::string name,
807 Mat *Adj,
int verb) {
816 auto p_miit = prb.find(name);
817 if (p_miit == prb.end()) {
819 "problem < %s > is not found (top tip: check spelling)",
823 std::vector<int> i_vec, j_vec;
824 j_vec.reserve(10000);
828 CHKERR PetscMalloc(i_vec.size() *
sizeof(
int), &_i);
829 CHKERR PetscMalloc(j_vec.size() *
sizeof(
int), &_j);
830 copy(i_vec.begin(), i_vec.end(), _i);
831 copy(j_vec.begin(), j_vec.end(), _j);
833 int nb_col_dofs = p_miit->getNbDofsCol();
834 CHKERR MatCreateMPIAdj(m_field.
get_comm(), i_vec.size() - 1, nb_col_dofs, _i,
835 _j, PETSC_NULL, Adj);
836 CHKERR MatSetOption(*Adj, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
844 const std::string name, Mat *Aij,
int verb) {
853 auto p_miit = prb.find(name);
854 if (p_miit == prb.end()) {
856 "problem < %s > is not found (top tip: check spelling)",
860 std::vector<int> i_vec, j_vec;
861 j_vec.reserve(10000);
865 int nb_local_dofs_row = p_miit->getNbLocalDofsRow();
866 int nb_local_dofs_col = p_miit->getNbLocalDofsCol();
869 CHKERR PetscMalloc(j_vec.size() *
sizeof(
double), &_a);
872 CHKERR ::MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, nb_local_dofs_row,
873 nb_local_dofs_col, &*i_vec.begin(),
874 &*j_vec.begin(), _a, &tmpMat);
875 CHKERR MatDuplicate(tmpMat, MAT_SHARE_NONZERO_PATTERN, Aij);
876 CHKERR MatDestroy(&tmpMat);
885 int row_print,
int col_print,
892 struct TestMatrixFillIn :
public FEMethod {
897 int rowPrint, colPrint;
899 TestMatrixFillIn(
CoreInterface *m_field_ptr, Mat
a,
int row_print,
901 : mFieldPtr(m_field_ptr),
A(
a), rowPrint(row_print),
902 colPrint(col_print){};
912 if (refinedFiniteElementsPtr->find(
913 numeredEntFiniteElementPtr->getEnt()) ==
914 refinedFiniteElementsPtr->end()) {
916 "data inconsistency");
919 auto row_dofs = getRowDofsPtr();
920 auto col_dofs = getColDofsPtr();
922 for (
auto cit = col_dofs->begin(); cit != col_dofs->end(); cit++) {
924 if (refinedEntitiesPtr->find((*cit)->getEnt()) ==
925 refinedEntitiesPtr->end()) {
927 "data inconsistency");
929 if (!(*cit)->getActive()) {
931 "data inconsistency");
934 FieldEntityEntFiniteElementAdjacencyMap_multiIndex::index<
937 boost::make_tuple((*cit)->getFieldEntityPtr()->getLocalUniqueId(),
938 numeredEntFiniteElementPtr->getLocalUniqueId()));
939 if (ait == adjacenciesPtr->end()) {
941 "adjacencies data inconsistency");
943 UId uid = ait->getEntUniqueId();
944 if (entitiesPtr->find(uid) == entitiesPtr->end()) {
946 "data inconsistency");
948 if (dofsPtr->find((*cit)->getLocalUniqueId()) == dofsPtr->end()) {
950 "data inconsistency");
954 if ((*cit)->getEntType() != MBVERTEX) {
957 col_dofs->get<
Ent_mi_tag>().equal_range((*cit)->getEnt());
958 int nb_dofs_on_ent = std::distance(range.first, range.second);
960 int max_order = (*cit)->getMaxOrder();
961 if ((*cit)->getNbOfCoeffs() * (*cit)->getOrderNbDofs(max_order) !=
969 <<
"Warning: Number of Dofs in Col different than number "
970 "of dofs for given entity order "
971 << (*cit)->getNbOfCoeffs() * (*cit)->getOrderNbDofs(max_order)
972 <<
" " << nb_dofs_on_ent;
977 for (
auto rit = row_dofs->begin(); rit != row_dofs->end(); rit++) {
979 if (refinedEntitiesPtr->find((*rit)->getEnt()) ==
980 refinedEntitiesPtr->end()) {
982 "data inconsistency");
984 if (!(*rit)->getActive()) {
986 "data inconsistency");
989 FieldEntityEntFiniteElementAdjacencyMap_multiIndex::index<
992 boost::make_tuple((*rit)->getFieldEntityPtr()->getLocalUniqueId(),
993 numeredEntFiniteElementPtr->getLocalUniqueId()));
994 if (ait == adjacenciesPtr->end()) {
996 MOFEM_LOG(
"SELF", Sev::error) << *(*rit);
997 MOFEM_LOG(
"SELF", Sev::error) << *(*rit);
998 MOFEM_LOG(
"SELF", Sev::error) << *numeredEntFiniteElementPtr;
999 MOFEM_LOG(
"SELF", Sev::error) <<
"dof: " << (*rit)->getBitRefLevel();
1001 <<
"fe: " << numeredEntFiniteElementPtr->getBitRefLevel();
1003 <<
"problem: " << problemPtr->getBitRefLevel();
1005 <<
"problem mask: " << problemPtr->getBitRefLevelMask();
1007 "adjacencies data inconsistency");
1009 UId uid = ait->getEntUniqueId();
1010 if (entitiesPtr->find(uid) == entitiesPtr->end()) {
1012 "data inconsistency");
1014 if (dofsPtr->find((*rit)->getLocalUniqueId()) == dofsPtr->end()) {
1016 "data inconsistency");
1019 int row = (*rit)->getPetscGlobalDofIdx();
1021 auto col_dofs = getColDofsPtr();
1022 for (
auto cit = col_dofs->begin(); cit != col_dofs->end(); cit++) {
1024 int col = (*cit)->getPetscGlobalDofIdx();
1026 if (row == rowPrint && col == colPrint) {
1027 MOFEM_LOG(
"SELF", Sev::noisy) <<
"fe:\n"
1028 << *numeredEntFiniteElementPtr;
1029 MOFEM_LOG(
"SELF", Sev::noisy) <<
"row:\n" << *(*rit);
1030 MOFEM_LOG(
"SELF", Sev::noisy) <<
"col:\n" << *(*cit);
1033 << numeredEntFiniteElementPtr->getBitRefLevel();
1034 MOFEM_LOG(
"SELF", Sev::noisy) <<
"row:\n"
1035 << (*rit)->getBitRefLevel();
1036 MOFEM_LOG(
"SELF", Sev::noisy) <<
"col:\n"
1037 << (*cit)->getBitRefLevel();
1040 CHKERR MatSetValue(
A, row, col, 1, INSERT_VALUES);
1043 if ((*rit)->getEntType() != MBVERTEX) {
1046 row_dofs->get<
Ent_mi_tag>().equal_range((*rit)->getEnt());
1047 int nb_dofs_on_ent = std::distance(range.first, range.second);
1049 int max_order = (*rit)->getMaxOrder();
1050 if ((*rit)->getNbOfCoeffs() * (*rit)->getOrderNbDofs(max_order) !=
1058 <<
"Warning: Number of Dofs in Row different than number "
1059 "of dofs for given entity order "
1060 << (*rit)->getNbOfCoeffs() * (*rit)->getOrderNbDofs(max_order)
1061 <<
" " << nb_dofs_on_ent;
1073 CHKERR MatAssemblyBegin(
A, MAT_FLUSH_ASSEMBLY);
1074 CHKERR MatAssemblyEnd(
A, MAT_FLUSH_ASSEMBLY);
1081 CHKERR MatSetOption(
A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
1084 MatView(
A, PETSC_VIEWER_STDOUT_WORLD);
1087 if (verb >=
NOISY) {
1088 MatView(
A, PETSC_VIEWER_DRAW_WORLD);
1093 TestMatrixFillIn method(&m_field,
A, row_print, col_print);
1098 auto p_miit = prb_set.find(problem_name);
1099 if (p_miit == prb_set.end())
1101 "problem < %s > not found (top tip: check spelling)",
1102 problem_name.c_str());
1103 MOFEM_LOG_C(
"WORLD", Sev::inform,
"check problem < %s >",
1104 problem_name.c_str());
1108 for (
auto &fe : *fe_ptr) {
1109 MOFEM_LOG_C(
"WORLD", Sev::verbose,
"\tcheck element %s",
1110 fe->getName().c_str());
1116 CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
1117 CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
1126MatrixManager::checkMPIAIJWithArraysMatrixFillIn<PetscGlobalIdx_mi_tag>(
1127 const std::string problem_name,
int row_print,
int col_print,
int verb) {
1131 CHKERR createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(problem_name,
A, verb);
1132 CHKERR MatSetOption(
A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
1139 const std::string problem_name,
int row_print,
int col_print,
int verb) {
1143 CHKERR createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_name,
A, verb);
1144 CHKERR MatSetOption(
A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
#define MatrixManagerFunctionBegin
Create adjacent matrices using different indices.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
#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 ...
@ 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 ...
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
Problem::EmptyFieldBlocks EmptyFieldBlocks
boost::shared_ptr< CacheTuple > CacheTupleSharedPtr
virtual MPI_Comm & get_comm() const =0
MPI_Comm & get_comm() const
PetscLogEvent MOFEM_EVENT_createMat
FieldEntity_multiIndex entsFields
entities on fields
int verbose
Verbosity level.
int sIze
MoFEM communicator size.
int rAnk
MOFEM communicator rank.
MPI_Comm mofemComm
MoFEM communicator.
std::reference_wrapper< moab::Interface > moab
moab database
FieldEntityEntFiniteElementAdjacencyMap_multiIndex entFEAdjacencies
adjacencies of elements to dofs
Create compressed matrix.
CreateRowComressedADJMatrix(moab::Interface &moab, MPI_Comm comm=PETSC_COMM_WORLD, int verbose=1)
MoFEMErrorCode getEntityAdjacenies(ProblemsByName::iterator p_miit, typename boost::multi_index::index< NumeredDofEntity_multiIndex, TAG >::type::iterator mit_row, boost::shared_ptr< FieldEntity > mofem_ent_ptr, std::vector< int > &dofs_col_view, int verb) const
Get element adjacencies.
FieldEntityEntFiniteElementAdjacencyMap_multiIndex::index< Unique_mi_tag >::type AdjByEnt
Problem_multiIndex::index< Problem_mi_tag >::type ProblemsByName
MoFEMErrorCode createMatArrays(ProblemsByName::iterator p_miit, const MatType type, std::vector< PetscInt > &i, std::vector< PetscInt > &j, const bool no_diagonals=true, int verb=QUIET) const
Create matrix adjacencies.
NumeredDofEntity_multiIndex::index< PetscGlobalIdx_mi_tag >::type DofByGlobalPetscIndex
structure for User Loop Methods on finite elements
Matrix manager is used to build and partition problems.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
PetscLogEvent MOFEM_EVENT_createSeqAIJCUSPARSEWithArrays
PetscLogEvent MOFEM_EVENT_createMPIAdjWithArrays
PetscLogEvent MOFEM_EVENT_checkMatrixFillIn
PetscLogEvent MOFEM_EVENT_createMPIAIJ
PetscLogEvent MOFEM_EVENT_createMPIAIJWithArrays
MoFEMErrorCode checkMatrixFillIn(const std::string problem_name, int row_print, int col_print, Mat A, int verb=QUIET)
check if matrix fill in correspond to finite element indices
PetscLogEvent MOFEM_EVENT_createSeqAIJWithArrays
MatrixManager(const MoFEM::Core &core)
PetscLogEvent MOFEM_EVENT_createMPIAIJCUSPARSEWithArrays
intrusive_ptr for managing petsc objects
base class for all interface classes