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<
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;
71 template <
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 auto check_block = [&empty_field_blocks](
auto &r_if,
auto &col_id) {
102 for(
auto &b : empty_field_blocks) {
103 if((b.first & r_if).any() && (b.second & col_id).any()) {
110 for (
auto &it :
r.first->entFePtr->getColFieldEnts()) {
111 if (
auto e = it.lock()) {
113 if (check_block((*mit_row)->getId(), e->getId())) {
115 if (
auto cache = e->entityCacheColDofs.lock()) {
116 const auto lo = cache->loHi[0];
117 const auto hi = cache->loHi[1];
118 for (
auto vit = lo; vit != hi; ++vit) {
120 const int idx = TAG::get_index(vit);
121 if (PetscLikely(idx >= 0))
122 dofs_col_view.push_back(idx);
125 const DofIdx nb_dofs_col = p_miit->getNbDofsCol();
126 if (PetscUnlikely(idx >= nb_dofs_col)) {
128 <<
"Problem with dof: " << std::endl
129 <<
"Rank " <<
rAnk <<
" : " << *(*vit);
132 "Index of dof larger than number of DOFs in column");
150 template <
typename TAG>
152 ProblemsByName::iterator p_miit,
const MatType
type,
153 std::vector<PetscInt> &
i, std::vector<PetscInt> &
j,
const bool no_diagonals,
159 auto cache = boost::make_shared<std::vector<EntityCacheNumeredDofs>>(
165 const auto uid = (*it)->getLocalUniqueId();
167 for (
auto lo =
r.first; lo !=
r.second; ++lo) {
169 if ((lo->getBitFEId() & p_miit->getBitFEId()).any()) {
171 const BitRefLevel &prb_bit = p_miit->getBitRefLevel();
172 const BitRefLevel &prb_mask = p_miit->getBitRefLevelMask();
173 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
176 if ((fe_bit & prb_mask) != fe_bit)
178 if ((fe_bit & prb_bit).none())
181 auto dit = p_miit->numeredColDofsPtr->lower_bound(uid);
183 decltype(dit) hi_dit;
184 if (dit != p_miit->numeredColDofsPtr->end())
185 hi_dit = p_miit->numeredColDofsPtr->upper_bound(
190 (*it)->entityCacheColDofs =
191 boost::shared_ptr<EntityCacheNumeredDofs>(cache, &((*cache)[idx]));
192 (*cache)[idx].loHi = {dit, hi_dit};
199 using NumeredDofEntitysByIdx =
204 const NumeredDofEntitysByIdx &dofs_row_by_idx =
205 p_miit->numeredRowDofsPtr->get<TAG>();
206 int nb_dofs_row = p_miit->getNbDofsRow();
207 if (nb_dofs_row == 0) {
209 p_miit->getName().c_str());
213 std::map<int, std::vector<int>> adjacent_dofs_on_other_parts;
220 TAG>::type::iterator miit_row,
223 if (TAG::IamNotPartitioned) {
228 CHKERR PetscLayoutSetBlockSize(layout, 1);
229 CHKERR PetscLayoutSetSize(layout, nb_dofs_row);
230 CHKERR PetscLayoutSetUp(layout);
231 PetscInt rstart, rend;
232 CHKERR PetscLayoutGetRange(layout, &rstart, &rend);
233 CHKERR PetscLayoutDestroy(&layout);
237 <<
"row lower " << rstart <<
" row upper " << rend;
241 miit_row = dofs_row_by_idx.lower_bound(rstart);
242 hi_miit_row = dofs_row_by_idx.lower_bound(rend);
243 if (std::distance(miit_row, hi_miit_row) != rend - rstart) {
246 "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
247 "rstart (%d != %d - %d = %d) ",
248 std::distance(miit_row, hi_miit_row), rend, rstart, rend - rstart);
258 std::vector<std::vector<int>> dofs_vec(
sIze);
260 boost::shared_ptr<FieldEntity> mofem_ent_ptr;
261 std::vector<int> dofs_col_view;
264 TAG>::type::iterator mit_row,
267 mit_row = dofs_row_by_idx.begin();
268 hi_mit_row = dofs_row_by_idx.end();
269 for (; mit_row != hi_mit_row; mit_row++) {
277 unsigned char pstatus = (*mit_row)->getPStatus();
278 if ((pstatus & PSTATUS_NOT_OWNED) &&
279 (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED))) {
281 bool get_adj_col =
true;
283 if (mofem_ent_ptr->getLocalUniqueId() ==
284 (*mit_row)->getFieldEntityPtr()->getLocalUniqueId()) {
291 mofem_ent_ptr = (*mit_row)->getFieldEntityPtr();
292 CHKERR getEntityAdjacenies<TAG>(p_miit, mit_row, mofem_ent_ptr,
293 dofs_col_view, verb);
296 std::sort(dofs_col_view.begin(), dofs_col_view.end());
297 std::vector<int>::iterator new_end =
298 std::unique(dofs_col_view.begin(), dofs_col_view.end());
299 int new_size = std::distance(dofs_col_view.begin(), new_end);
300 dofs_col_view.resize(new_size);
304 int owner = (*mit_row)->getOwnerProc();
305 dofs_vec[owner].emplace_back(TAG::get_index(mit_row));
306 dofs_vec[owner].emplace_back(
307 dofs_col_view.size());
309 dofs_vec[owner].insert(dofs_vec[owner].end(), dofs_col_view.begin(),
310 dofs_col_view.end());
316 std::vector<int> dofs_vec_length(
sIze);
317 for (
int proc = 0; proc <
sIze; proc++) {
319 if (!dofs_vec[proc].empty()) {
321 dofs_vec_length[proc] = dofs_vec[proc].size();
326 dofs_vec_length[proc] = 0;
330 std::vector<MPI_Status> status(
sIze);
334 CHKERR PetscGatherNumberOfMessages(comm, NULL, &dofs_vec_length[0],
341 CHKERR PetscGatherMessageLengths(comm, nsends, nrecvs, &dofs_vec_length[0],
346 CHKERR PetscCommGetNewTag(comm, &tag);
351 MPI_Request *r_waits;
356 CHKERR PetscPostIrecvInt(comm, tag, nrecvs, onodes, olengths, &rbuf,
359 MPI_Request *s_waits;
360 CHKERR PetscMalloc1(nsends, &s_waits);
363 for (
int proc = 0, kk = 0; proc <
sIze; proc++) {
364 if (!dofs_vec_length[proc])
366 CHKERR MPI_Isend(&(dofs_vec[proc])[0],
367 dofs_vec_length[proc],
369 tag, comm, s_waits + kk);
375 CHKERR MPI_Waitall(nrecvs, r_waits, &status[0]);
379 CHKERR MPI_Waitall(nsends, s_waits, &status[0]);
382 for (
int kk = 0; kk < nrecvs; kk++) {
384 int len = olengths[kk];
385 int *data_from_proc = rbuf[kk];
387 for (
int ii = 0; ii < len;) {
389 int row_idx = data_from_proc[ii++];
390 int nb_adj_dofs = data_from_proc[ii++];
394 DofByGlobalPetscIndex::iterator dit;
400 "dof %d can not be found in problem", row_idx);
404 for (
int jj = 0; jj < nb_adj_dofs; jj++) {
405 adjacent_dofs_on_other_parts[row_idx].push_back(data_from_proc[ii++]);
411 CHKERR PetscFree(s_waits);
412 CHKERR PetscFree(rbuf[0]);
414 CHKERR PetscFree(r_waits);
416 CHKERR PetscFree(olengths);
418 miit_row = dofs_row_by_idx.begin();
419 hi_miit_row = dofs_row_by_idx.end();
421 CHKERR PetscCommDestroy(&comm);
424 boost::shared_ptr<FieldEntity> mofem_ent_ptr;
425 int row_last_evaluated_idx = -1;
427 std::vector<int> dofs_vec;
428 std::vector<int> dofs_col_view;
431 int nb_loc_row_from_iterators = 0;
432 unsigned int rows_to_fill = p_miit->getNbLocalDofsRow();
433 i.reserve(rows_to_fill + 1);
434 for (; miit_row != hi_miit_row; miit_row++) {
436 if (!TAG::IamNotPartitioned) {
437 if (
static_cast<int>((*miit_row)->getPart()) !=
rAnk)
441 nb_loc_row_from_iterators++;
444 i.push_back(
j.size());
451 : (mofem_ent_ptr->getLocalUniqueId() !=
452 (*miit_row)->getFieldEntityPtr()->getLocalUniqueId())) {
455 mofem_ent_ptr = (*miit_row)->getFieldEntityPtr();
456 CHKERR getEntityAdjacenies<TAG>(p_miit, miit_row, mofem_ent_ptr,
457 dofs_col_view, verb);
458 row_last_evaluated_idx = TAG::get_index(miit_row);
462 dofs_vec.insert(dofs_vec.end(), dofs_col_view.begin(),
463 dofs_col_view.end());
465 unsigned char pstatus = (*miit_row)->getPStatus();
467 std::map<int, std::vector<int>>::iterator mit;
468 mit = adjacent_dofs_on_other_parts.find(row_last_evaluated_idx);
469 if (mit == adjacent_dofs_on_other_parts.end()) {
477 dofs_vec.insert(dofs_vec.end(), mit->second.begin(),
483 sort(dofs_vec.begin(), dofs_vec.end());
484 std::vector<int>::iterator new_end =
485 unique(dofs_vec.begin(), dofs_vec.end());
486 int new_size = std::distance(dofs_vec.begin(), new_end);
487 dofs_vec.resize(new_size);
491 if (
j.capacity() <
j.size() + dofs_vec.size()) {
494 unsigned int nb_nonzero =
j.size() + dofs_vec.size();
495 unsigned int average_row_fill = nb_nonzero /
i.size() + 1;
496 j.reserve(rows_to_fill * average_row_fill);
499 auto hi_diit = dofs_vec.end();
500 for (
auto diit = dofs_vec.begin(); diit != hi_diit; diit++) {
503 if (*diit == TAG::get_index(miit_row)) {
512 i.push_back(
j.size());
514 if (strcmp(
type, MATMPIADJ) == 0) {
517 if (
i.size() - 1 != (
unsigned int)nb_loc_row_from_iterators) {
520 "Number of rows from iterator is different than size of rows in "
522 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
524 (
unsigned int)nb_loc_row_from_iterators,
i.size() - 1);
527 }
else if (strcmp(
type, MATMPIAIJ) == 0) {
530 if (
i.size() - 1 != (
unsigned int)nb_loc_row_from_iterators) {
533 "Number of rows from iterator is different than size of rows in "
535 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
537 (
unsigned int)nb_loc_row_from_iterators,
i.size() - 1);
539 PetscInt nb_local_dofs_row = p_miit->getNbLocalDofsRow();
540 if ((
unsigned int)nb_local_dofs_row !=
i.size() - 1) {
543 "Number of rows is different than size of rows in compressed row "
544 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
546 (
unsigned int)nb_local_dofs_row,
i.size() - 1);
549 }
else if (strcmp(
type, MATAIJ) == 0) {
552 if (
i.size() - 1 != (
unsigned int)nb_loc_row_from_iterators) {
555 "Number of rows form iterator is different than size of rows in "
557 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
559 (
unsigned int)nb_loc_row_from_iterators,
i.size() - 1);
561 PetscInt nb_local_dofs_row = p_miit->getNbLocalDofsRow();
562 if ((
unsigned int)nb_local_dofs_row !=
i.size() - 1) {
565 "Number of rows is different than size of rows in compressed row "
566 "matrix (unsigned int)nb_local_dofs_row != i.size() - 1, i.e. %d != "
568 (
unsigned int)nb_local_dofs_row,
i.size() - 1);
571 }
else if (strcmp(
type, MATAIJCUSPARSE) == 0) {
573 if (
i.size() - 1 != (
unsigned int)nb_loc_row_from_iterators) {
574 SETERRQ(
get_comm(), PETSC_ERR_ARG_SIZ,
"data inconsistency");
576 PetscInt nb_local_dofs_row = p_miit->getNbLocalDofsRow();
577 if ((
unsigned int)nb_local_dofs_row !=
i.size() - 1) {
578 SETERRQ(
get_comm(), PETSC_ERR_ARG_SIZ,
"data inconsistency");
583 SETERRQ(
get_comm(), PETSC_ERR_ARG_NULL,
"not implemented");
598 : cOre(const_cast<
MoFEM::
Core &>(core)) {
599 PetscLogEventRegister(
"MatrixManagerCreateMPIAIJ", 0,
601 PetscLogEventRegister(
"MatrixManagerCreateMPIAIJWithArrays", 0,
603 PetscLogEventRegister(
"MatrixManagerCreateMPIAdjWithArrays", 0,
605 PetscLogEventRegister(
"MatrixManagerCreateMPIAIJCUSPARSEWithArrays", 0,
607 PetscLogEventRegister(
"MatrixManagerCreateSeqAIJCUSPARSEWithArrays", 0,
609 PetscLogEventRegister(
"MatrixManagerCreateSeqAIJWithArrays", 0,
611 PetscLogEventRegister(
"MatrixManagerCheckMPIAIJWithArraysMatrixFillIn", 0,
617 const std::string name, Mat *Aij,
int verb) {
627 auto p_miit = prb.find(name);
628 if (p_miit == prb.end()) {
630 "problem < %s > is not found (top tip: check spelling)",
634 std::vector<int> i_vec, j_vec;
635 j_vec.reserve(10000);
637 p_miit, MATMPIAIJ, i_vec, j_vec,
false, verb);
639 int nb_row_dofs = p_miit->getNbDofsRow();
640 int nb_col_dofs = p_miit->getNbDofsCol();
641 int nb_local_dofs_row = p_miit->getNbLocalDofsRow();
642 int nb_local_dofs_col = p_miit->getNbLocalDofsCol();
644 CHKERR ::MatCreateMPIAIJWithArrays(
645 m_field.
get_comm(), nb_local_dofs_row, nb_local_dofs_col, nb_row_dofs,
646 nb_col_dofs, &*i_vec.begin(), &*j_vec.begin(), PETSC_NULL, Aij);
654 MatrixManager::createMPIAIJCUSPARSEWithArrays<PetscGlobalIdx_mi_tag>(
655 const std::string name, Mat *Aij,
int verb) {
665 auto p_miit = prb.find(name);
666 if (p_miit == prb.end()) {
668 "problem < %s > is not found (top tip: check spelling)",
672 std::vector<int> i_vec, j_vec;
673 j_vec.reserve(10000);
675 p_miit, MATAIJCUSPARSE, i_vec, j_vec,
false, verb);
677 int nb_local_dofs_row = p_miit->getNbLocalDofsRow();
678 int nb_local_dofs_col = p_miit->getNbLocalDofsCol();
680 auto get_layout = [&]() {
681 int start_ranges, end_ranges;
684 CHKERR PetscLayoutSetBlockSize(layout, 1);
685 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs_col);
686 CHKERR PetscLayoutSetUp(layout);
687 CHKERR PetscLayoutGetRange(layout, &start_ranges, &end_ranges);
688 CHKERR PetscLayoutDestroy(&layout);
689 return std::make_pair(start_ranges, end_ranges);
692 auto get_nnz = [&](
auto &d_nnz,
auto &o_nnz) {
694 auto layout = get_layout();
696 for (
int i = 0;
i != nb_local_dofs_row; ++
i) {
697 for (;
j != i_vec[
i + 1]; ++
j) {
698 if (j_vec[
j] < layout.second && j_vec[
j] >= layout.first)
707 std::vector<int> d_nnz(nb_local_dofs_row, 0), o_nnz(nb_local_dofs_row, 0);
708 CHKERR get_nnz(d_nnz, o_nnz);
710 #ifdef PETSC_HAVE_CUDA
711 CHKERR ::MatCreateAIJCUSPARSE(m_field.
get_comm(), nb_local_dofs_row,
712 nb_local_dofs_col, nb_row_dofs, nb_col_dofs, 0,
713 &*d_nnz.begin(), 0, &*o_nnz.begin(), Aij);
716 "Error: To use this matrix type compile PETSc with CUDA.");
726 MatrixManager::createSeqAIJCUSPARSEWithArrays<PetscLocalIdx_mi_tag>(
727 const std::string name, Mat *Aij,
int verb) {
730 #ifdef PETSC_HAVE_CUDA
734 "Not implemented type of matrix yet, try MPI version (aijcusparse)");
742 MatrixManager::createMPIAIJ<PetscGlobalIdx_mi_tag>(
const std::string name,
743 Mat *Aij,
int verb) {
752 auto p_miit = prb.find(name);
753 if (p_miit == prb.end()) {
755 "problem < %s > is not found (top tip: check spelling)",
759 std::vector<int> i_vec, j_vec;
760 j_vec.reserve(10000);
762 p_miit, MATMPIAIJ, i_vec, j_vec,
false, verb);
764 int nb_row_dofs = p_miit->getNbDofsRow();
765 int nb_col_dofs = p_miit->getNbDofsCol();
766 int nb_local_dofs_row = p_miit->getNbLocalDofsRow();
767 int nb_local_dofs_col = p_miit->getNbLocalDofsCol();
769 auto get_layout = [&]() {
770 int start_ranges, end_ranges;
773 CHKERR PetscLayoutSetBlockSize(layout, 1);
774 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs_col);
775 CHKERR PetscLayoutSetUp(layout);
776 CHKERR PetscLayoutGetRange(layout, &start_ranges, &end_ranges);
777 CHKERR PetscLayoutDestroy(&layout);
778 return std::make_pair(start_ranges, end_ranges);
781 auto get_nnz = [&](
auto &d_nnz,
auto &o_nnz) {
783 auto layout = get_layout();
785 for (
int i = 0;
i != nb_local_dofs_row; ++
i) {
786 for (;
j != i_vec[
i + 1]; ++
j) {
787 if (j_vec[
j] < layout.second && j_vec[
j] >= layout.first)
796 std::vector<int> d_nnz(nb_local_dofs_row, 0), o_nnz(nb_local_dofs_row, 0);
797 CHKERR get_nnz(d_nnz, o_nnz);
800 CHKERR MatSetSizes(*Aij, nb_local_dofs_row, nb_local_dofs_col, nb_row_dofs,
802 CHKERR MatSetType(*Aij, MATMPIAIJ);
803 CHKERR MatMPIAIJSetPreallocation(*Aij, 0, &*d_nnz.begin(), 0,
812 MatrixManager::createMPIAdjWithArrays<Idx_mi_tag>(
const std::string name,
813 Mat *Adj,
int verb) {
822 auto p_miit = prb.find(name);
823 if (p_miit == prb.end()) {
825 "problem < %s > is not found (top tip: check spelling)",
829 std::vector<int> i_vec, j_vec;
830 j_vec.reserve(10000);
834 CHKERR PetscMalloc(i_vec.size() *
sizeof(
int), &_i);
835 CHKERR PetscMalloc(j_vec.size() *
sizeof(
int), &_j);
836 copy(i_vec.begin(), i_vec.end(), _i);
837 copy(j_vec.begin(), j_vec.end(), _j);
839 int nb_col_dofs = p_miit->getNbDofsCol();
840 CHKERR MatCreateMPIAdj(m_field.
get_comm(), i_vec.size() - 1, nb_col_dofs, _i,
841 _j, PETSC_NULL, Adj);
842 CHKERR MatSetOption(*Adj, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
850 const std::string name, Mat *Aij,
int verb) {
859 auto p_miit = prb.find(name);
860 if (p_miit == prb.end()) {
862 "problem < %s > is not found (top tip: check spelling)",
866 std::vector<int> i_vec, j_vec;
867 j_vec.reserve(10000);
871 int nb_local_dofs_row = p_miit->getNbLocalDofsRow();
872 int nb_local_dofs_col = p_miit->getNbLocalDofsCol();
875 CHKERR PetscMalloc(j_vec.size() *
sizeof(
double), &_a);
878 CHKERR ::MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, nb_local_dofs_row,
879 nb_local_dofs_col, &*i_vec.begin(),
880 &*j_vec.begin(), _a, &tmpMat);
881 CHKERR MatDuplicate(tmpMat, MAT_SHARE_NONZERO_PATTERN, Aij);
882 CHKERR MatDestroy(&tmpMat);
891 int row_print,
int col_print,
898 struct TestMatrixFillIn :
public FEMethod {
903 int rowPrint, colPrint;
905 TestMatrixFillIn(
CoreInterface *m_field_ptr, Mat
a,
int row_print,
907 : mFieldPtr(m_field_ptr),
A(
a), rowPrint(row_print),
908 colPrint(col_print){};
918 if (refinedFiniteElementsPtr->find(
919 numeredEntFiniteElementPtr->getEnt()) ==
920 refinedFiniteElementsPtr->end()) {
922 "data inconsistency");
925 auto row_dofs = getRowDofsPtr();
926 auto col_dofs = getColDofsPtr();
928 for (
auto cit = col_dofs->begin(); cit != col_dofs->end(); cit++) {
930 if (refinedEntitiesPtr->find((*cit)->getEnt()) ==
931 refinedEntitiesPtr->end()) {
933 "data inconsistency");
935 if (!(*cit)->getActive()) {
937 "data inconsistency");
940 FieldEntityEntFiniteElementAdjacencyMap_multiIndex::index<
943 boost::make_tuple((*cit)->getFieldEntityPtr()->getLocalUniqueId(),
944 numeredEntFiniteElementPtr->getLocalUniqueId()));
945 if (ait == adjacenciesPtr->end()) {
947 "adjacencies data inconsistency");
949 UId uid = ait->getEntUniqueId();
950 if (entitiesPtr->find(uid) == entitiesPtr->end()) {
952 "data inconsistency");
954 if (dofsPtr->find((*cit)->getLocalUniqueId()) == dofsPtr->end()) {
956 "data inconsistency");
960 if ((*cit)->getEntType() != MBVERTEX) {
963 col_dofs->get<
Ent_mi_tag>().equal_range((*cit)->getEnt());
964 int nb_dofs_on_ent = std::distance(range.first, range.second);
966 int max_order = (*cit)->getMaxOrder();
967 if ((*cit)->getNbOfCoeffs() * (*cit)->getOrderNbDofs(max_order) !=
975 <<
"Warning: Number of Dofs in Col different than number "
976 "of dofs for given entity order "
977 << (*cit)->getNbOfCoeffs() * (*cit)->getOrderNbDofs(max_order)
978 <<
" " << nb_dofs_on_ent;
983 for (
auto rit = row_dofs->begin(); rit != row_dofs->end(); rit++) {
985 if (refinedEntitiesPtr->find((*rit)->getEnt()) ==
986 refinedEntitiesPtr->end()) {
988 "data inconsistency");
990 if (!(*rit)->getActive()) {
992 "data inconsistency");
995 FieldEntityEntFiniteElementAdjacencyMap_multiIndex::index<
998 boost::make_tuple((*rit)->getFieldEntityPtr()->getLocalUniqueId(),
999 numeredEntFiniteElementPtr->getLocalUniqueId()));
1000 if (ait == adjacenciesPtr->end()) {
1002 MOFEM_LOG(
"SELF", Sev::error) << *(*rit);
1003 MOFEM_LOG(
"SELF", Sev::error) << *(*rit);
1004 MOFEM_LOG(
"SELF", Sev::error) << *numeredEntFiniteElementPtr;
1005 MOFEM_LOG(
"SELF", Sev::error) <<
"dof: " << (*rit)->getBitRefLevel();
1007 <<
"fe: " << numeredEntFiniteElementPtr->getBitRefLevel();
1009 <<
"problem: " << problemPtr->getBitRefLevel();
1011 <<
"problem mask: " << problemPtr->getBitRefLevelMask();
1013 "adjacencies data inconsistency");
1015 UId uid = ait->getEntUniqueId();
1016 if (entitiesPtr->find(uid) == entitiesPtr->end()) {
1018 "data inconsistency");
1020 if (dofsPtr->find((*rit)->getLocalUniqueId()) == dofsPtr->end()) {
1022 "data inconsistency");
1025 int row = (*rit)->getPetscGlobalDofIdx();
1027 auto col_dofs = getColDofsPtr();
1028 for (
auto cit = col_dofs->begin(); cit != col_dofs->end(); cit++) {
1030 int col = (*cit)->getPetscGlobalDofIdx();
1032 if (row == rowPrint && col == colPrint) {
1033 MOFEM_LOG(
"SELF", Sev::noisy) <<
"fe:\n"
1034 << *numeredEntFiniteElementPtr;
1035 MOFEM_LOG(
"SELF", Sev::noisy) <<
"row:\n" << *(*rit);
1036 MOFEM_LOG(
"SELF", Sev::noisy) <<
"col:\n" << *(*cit);
1039 << numeredEntFiniteElementPtr->getBitRefLevel();
1040 MOFEM_LOG(
"SELF", Sev::noisy) <<
"row:\n"
1041 << (*rit)->getBitRefLevel();
1042 MOFEM_LOG(
"SELF", Sev::noisy) <<
"col:\n"
1043 << (*cit)->getBitRefLevel();
1046 CHKERR MatSetValue(
A, row, col, 1, INSERT_VALUES);
1049 if ((*rit)->getEntType() != MBVERTEX) {
1052 row_dofs->get<
Ent_mi_tag>().equal_range((*rit)->getEnt());
1053 int nb_dofs_on_ent = std::distance(range.first, range.second);
1055 int max_order = (*rit)->getMaxOrder();
1056 if ((*rit)->getNbOfCoeffs() * (*rit)->getOrderNbDofs(max_order) !=
1064 <<
"Warning: Number of Dofs in Row different than number "
1065 "of dofs for given entity order "
1066 << (*rit)->getNbOfCoeffs() * (*rit)->getOrderNbDofs(max_order)
1067 <<
" " << nb_dofs_on_ent;
1079 CHKERR MatAssemblyBegin(
A, MAT_FLUSH_ASSEMBLY);
1080 CHKERR MatAssemblyEnd(
A, MAT_FLUSH_ASSEMBLY);
1087 CHKERR MatSetOption(
A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
1090 MatView(
A, PETSC_VIEWER_STDOUT_WORLD);
1093 if (verb >=
NOISY) {
1094 MatView(
A, PETSC_VIEWER_DRAW_WORLD);
1099 TestMatrixFillIn method(&m_field,
A, row_print, col_print);
1104 auto p_miit = prb_set.find(problem_name);
1105 if (p_miit == prb_set.end())
1107 "problem < %s > not found (top tip: check spelling)",
1108 problem_name.c_str());
1109 MOFEM_LOG_C(
"WORLD", Sev::inform,
"check problem < %s >",
1110 problem_name.c_str());
1114 for (
auto &fe : *fe_ptr) {
1115 MOFEM_LOG_C(
"WORLD", Sev::verbose,
"\tcheck element %s",
1116 fe->getName().c_str());
1122 CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
1123 CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
1132 MatrixManager::checkMPIAIJWithArraysMatrixFillIn<PetscGlobalIdx_mi_tag>(
1133 const std::string problem_name,
int row_print,
int col_print,
int verb) {
1137 CHKERR createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(problem_name,
A, verb);
1138 CHKERR MatSetOption(
A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
1145 const std::string problem_name,
int row_print,
int col_print,
int verb) {
1149 CHKERR createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_name,
A, verb);
1150 CHKERR MatSetOption(
A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);