v0.14.0
Loading...
Searching...
No Matches
ProblemsManager

Adding and managing problems. More...

Collaboration diagram for ProblemsManager:

Files

file  ProblemsManager.cpp
 Managing complexities for problem.
 
file  ProblemsManager.hpp
 Interface managing problems.
 

Classes

struct  MoFEM::ProblemsManager
 Problem manager is used to build and partition problems. More...
 

Functions

DEPRECATED MoFEMErrorCode MoFEM::ProblemsManager::partitionMesh (const Range &ents, const int dim, const int adj_dim, const int n_parts, Tag *th_vertex_weights=nullptr, Tag *th_edge_weights=nullptr, Tag *th_part_weights=nullptr, int verb=VERBOSE, const bool debug=false)
 Set partition tag to each finite element in the problem.
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblem (const std::string name, const bool square_matrix, int verb=VERBOSE)
 build problem data structures
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblem (Problem *problem_ptr, const bool square_matrix, int verb=VERBOSE)
 build problem data structures
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh (const std::string name, const bool square_matrix, int verb=VERBOSE)
 build problem data structures, assuming that mesh is distributed (collective)
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh (Problem *problem_ptr, const bool square_matrix=true, int verb=VERBOSE)
 build problem data structures, assuming that mesh is distributed (collective)
 
MoFEMErrorCode MoFEM::ProblemsManager::buildSubProblem (const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, boost::shared_ptr< Range > > *entityMapRow=nullptr, const map< std::string, boost::shared_ptr< Range > > *entityMapCol=nullptr, int verb=VERBOSE)
 build sub problem
 
MoFEMErrorCode MoFEM::ProblemsManager::buildComposedProblem (const std::string out_name, const std::vector< std::string > add_row_problems, const std::vector< std::string > add_col_problems, const bool square_matrix=true, int verb=1)
 build composite problem
 
MoFEMErrorCode MoFEM::ProblemsManager::inheritPartition (const std::string name, const std::string problem_for_rows, bool copy_rows, const std::string problem_for_cols, bool copy_cols, int verb=VERBOSE)
 build indexing and partition problem inheriting indexing and partitioning from two other problems
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionSimpleProblem (const std::string name, int verb=VERBOSE)
 partition problem dofs
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionProblem (const std::string name, int verb=VERBOSE)
 partition problem dofs (collective)
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionFiniteElements (const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
 partition finite elements
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofs (const std::string name, int verb=VERBOSE)
 determine ghost nodes
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofsOnDistributedMesh (const std::string name, int verb=VERBOSE)
 determine ghost nodes on distributed meshes
 
MoFEMErrorCode MoFEM::ProblemsManager::getFEMeshset (const std::string prb_name, const std::string &fe_name, EntityHandle *meshset) const
 create add entities of finite element in the problem
 
MoFEMErrorCode MoFEM::ProblemsManager::getProblemElementsLayout (const std::string name, const std::string &fe_name, PetscLayout *layout) const
 Get layout of elements in the problem.
 
MoFEMErrorCode MoFEM::ProblemsManager::removeDofsOnEntities (const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
 Remove DOFs from problem.
 
MoFEMErrorCode MoFEM::ProblemsManager::removeDofsOnEntities (const std::string problem_name, const std::string field_name, const BitRefLevel bit_ref_level, const BitRefLevel bit_ref_mask, Range *ents_ptr=nullptr, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
 Remove DOFs from problem by bit ref level.
 

Read load and boradcoast

@

MoFEMErrorCode MoFEM::CommInterface::partitionMesh (const Range &ents, const int dim, const int adj_dim, const int n_parts, Tag *th_vertex_weights=nullptr, Tag *th_edge_weights=nullptr, Tag *th_part_weights=nullptr, int verb=VERBOSE, const bool debug=false)
 Set partition tag to each finite element in the problem.
 

Detailed Description

Adding and managing problems.

Function Documentation

◆ buildComposedProblem()

MoFEMErrorCode MoFEM::ProblemsManager::buildComposedProblem ( const std::string out_name,
const std::vector< std::string > add_row_problems,
const std::vector< std::string > add_col_problems,
const bool square_matrix = true,
int verb = 1 )

build composite problem

Parameters
out_namename of build problem
add_row_problemsvector of add row problems
add_col_problemsvector of add col problems
square_matrixtrue if structurally squared matrix
verbverbosity level
Returns
error code

Definition at line 1250 of file ProblemsManager.cpp.

1253 {
1255 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1257 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1259 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1260 MoFEM::Interface &m_field = cOre;
1261 auto problems_ptr = m_field.get_problems();
1263
1264 CHKERR m_field.clear_problem(out_name);
1265 // get reference to all problems
1266 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1267 ProblemByName &problems_by_name =
1268 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1269
1270 // Get iterators to out problem, i.e. build problem
1271 ProblemByName::iterator out_problem_it = problems_by_name.find(out_name);
1272 if (out_problem_it == problems_by_name.end()) {
1273 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1274 "problem with name < %s > not defined (top tip check spelling)",
1275 out_name.c_str());
1276 }
1277 // Make data structure for composed-problem data
1278 out_problem_it->composedProblemsData =
1279 boost::make_shared<ComposedProblemsData>();
1280 boost::shared_ptr<ComposedProblemsData> cmp_prb_data =
1281 out_problem_it->getComposedProblemsData();
1282
1283 const std::vector<std::string> *add_prb[] = {&add_row_problems,
1284 &add_col_problems};
1285 std::vector<const Problem *> *add_prb_ptr[] = {&cmp_prb_data->rowProblemsAdd,
1286 &cmp_prb_data->colProblemsAdd};
1287 std::vector<SmartPetscObj<IS>> *add_prb_is[] = {&cmp_prb_data->rowIs,
1288 &cmp_prb_data->colIs};
1289
1290 // Get local indices counter
1291 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1292 &out_problem_it->nbLocDofsCol};
1293 // Get global indices counter
1294 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1295
1296 // Set number of ghost nodes to zero
1297 {
1298 out_problem_it->nbDofsRow = 0;
1299 out_problem_it->nbDofsCol = 0;
1300 out_problem_it->nbLocDofsRow = 0;
1301 out_problem_it->nbLocDofsCol = 0;
1302 out_problem_it->nbGhostDofsRow = 0;
1303 out_problem_it->nbGhostDofsCol = 0;
1304 }
1305 int nb_dofs_reserve[] = {0, 0};
1306
1307 // Loop over rows and columns in the main problem and sub-problems
1308 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1309 add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1310 add_prb_is[ss]->reserve(add_prb[ss]->size());
1311 for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1312 vit != add_prb[ss]->end(); vit++) {
1313 ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1314 if (prb_it == problems_by_name.end()) {
1315 SETERRQ1(
1316 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1317 "problem with name < %s > not defined (top tip check spelling)",
1318 vit->c_str());
1319 }
1320 add_prb_ptr[ss]->push_back(&*prb_it);
1321 // set number of dofs on rows and columns
1322 if (ss == 0) {
1323 // row
1324 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1325 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1326 nb_dofs_reserve[ss] +=
1327 add_prb_ptr[ss]->back()->numeredRowDofsPtr->size();
1328 } else {
1329 // column
1330 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1331 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1332 nb_dofs_reserve[ss] +=
1333 add_prb_ptr[ss]->back()->numeredColDofsPtr->size();
1334 }
1335 }
1336 }
1337 // if squre problem, rows and columns are the same
1338 if (square_matrix) {
1339 add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1340 add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1341 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1342 *nb_dofs[1] = *nb_dofs[0];
1343 *nb_local_dofs[1] = *nb_local_dofs[0];
1344 }
1345
1346 // reserve memory for dofs
1347 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1348 // Reserve memory
1349 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1350 dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1351 dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1352 if (!ss)
1353 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1354 else
1355 out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1356 }
1357
1358 // Push back DOFs
1359 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1360 NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1361 dit,
1362 hi_dit;
1363 int shift_glob = 0;
1364 int shift_loc = 0;
1365 for (unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1366 PetscInt *dofs_out_idx_ptr;
1367 int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1368 CHKERR PetscMalloc(nb_local_dofs * sizeof(int), &dofs_out_idx_ptr);
1369 if (ss == 0) {
1370 dit = (*add_prb_ptr[ss])[pp]
1371 ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1372 .begin();
1373 hi_dit = (*add_prb_ptr[ss])[pp]
1374 ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1375 .end();
1376 } else {
1377 dit = (*add_prb_ptr[ss])[pp]
1378 ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1379 .begin();
1380 hi_dit = (*add_prb_ptr[ss])[pp]
1381 ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1382 .end();
1383 }
1384 int is_nb = 0;
1385 for (; dit != hi_dit; dit++) {
1386 const BitRefLevel &prb_bit = out_problem_it->getBitRefLevel();
1387 const BitRefLevel &prb_mask = out_problem_it->getBitRefLevelMask();
1388 const BitRefLevel &dof_bit = dit->get()->getBitRefLevel();
1389 if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1390 continue;
1391 const int rank = m_field.get_comm_rank();
1392 const int part = dit->get()->getPart();
1393 const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1394 const int loc_idx =
1395 (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1396 : -1;
1397 dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1398 glob_idx, loc_idx, part);
1399 if (part == rank) {
1400 dofs_out_idx_ptr[is_nb++] = glob_idx;
1401 }
1402 }
1403 if (is_nb > nb_local_dofs) {
1404 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1405 "Data inconsistency");
1406 }
1407 IS is;
1408 CHKERR ISCreateGeneral(m_field.get_comm(), is_nb, dofs_out_idx_ptr,
1409 PETSC_OWN_POINTER, &is);
1410 auto smart_is = SmartPetscObj<IS>(is);
1411 (*add_prb_is[ss]).push_back(smart_is);
1412 if (ss == 0) {
1413 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1414 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1415 } else {
1416 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1417 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1418 }
1419 if (square_matrix) {
1420 (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1421 (*add_prb_is[1]).push_back(smart_is);
1422 }
1423 }
1424 }
1425
1426 if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1427 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1428 }
1429
1430 // Insert DOFs to problem multi-index
1431 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1432 auto hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->end()
1433 : out_problem_it->numeredColDofsPtr->end();
1434 for (auto &v : *dofs_array[ss])
1435 hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->emplace_hint(
1436 hint, dofs_array[ss], &v)
1437 : out_problem_it->numeredColDofsPtr->emplace_hint(
1438 hint, dofs_array[ss], &v);
1439 }
1440
1441 // Compress DOFs
1442 *nb_dofs[0] = 0;
1443 *nb_dofs[1] = 0;
1444 *nb_local_dofs[0] = 0;
1445 *nb_local_dofs[1] = 0;
1446 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1447
1448 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1449 if (ss == 0) {
1450 dofs_ptr = out_problem_it->numeredRowDofsPtr;
1451 } else {
1452 dofs_ptr = out_problem_it->numeredColDofsPtr;
1453 }
1454 NumeredDofEntityByUId::iterator dit, hi_dit;
1455 dit = dofs_ptr->get<Unique_mi_tag>().begin();
1456 hi_dit = dofs_ptr->get<Unique_mi_tag>().end();
1457 std::vector<int> idx;
1458 idx.reserve(std::distance(dit, hi_dit));
1459 // set dofs in order entity and dof number on entity
1460 for (; dit != hi_dit; dit++) {
1461 if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1462 bool success = dofs_ptr->get<Unique_mi_tag>().modify(
1463 dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1464 if (!success) {
1465 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1466 "modification unsuccessful");
1467 }
1468 idx.push_back(dit->get()->getPetscGlobalDofIdx());
1469 } else {
1470 if (dit->get()->getPetscLocalDofIdx() != -1) {
1471 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1472 "local index should be negative");
1473 }
1474 }
1475 }
1476 if (square_matrix) {
1477 *nb_local_dofs[1] = *nb_local_dofs[0];
1478 }
1479
1480 // set new dofs mapping
1481 IS is;
1482 CHKERR ISCreateGeneral(m_field.get_comm(), idx.size(), &*idx.begin(),
1483 PETSC_USE_POINTER, &is);
1484 CHKERR ISGetSize(is, nb_dofs[ss]);
1485 if (square_matrix) {
1486 *nb_dofs[1] = *nb_dofs[0];
1487 }
1488
1489 AO ao;
1490 CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1491 for (unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1492 CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1493
1494 // Set DOFs numeration
1495 {
1496 std::vector<int> idx_new;
1497 idx_new.reserve(dofs_ptr->size());
1498 for (NumeredDofEntityByUId::iterator dit =
1499 dofs_ptr->get<Unique_mi_tag>().begin();
1500 dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1501 idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1502 }
1503 // set new global dofs numeration
1504 IS is_new;
1505 CHKERR ISCreateGeneral(m_field.get_comm(), idx_new.size(),
1506 &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1507 CHKERR AOApplicationToPetscIS(ao, is_new);
1508 // set global indices to multi-index
1509 std::vector<int>::iterator vit = idx_new.begin();
1510 for (NumeredDofEntityByUId::iterator dit =
1511 dofs_ptr->get<Unique_mi_tag>().begin();
1512 dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1513 bool success =
1514 dofs_ptr->modify(dit, NumeredDofEntity_part_and_glob_idx_change(
1515 dit->get()->getPart(), *(vit++)));
1516 if (!success) {
1517 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1518 "modification unsuccessful");
1519 }
1520 }
1521 CHKERR ISDestroy(&is_new);
1522 }
1523 CHKERR ISDestroy(&is);
1524 CHKERR AODestroy(&ao);
1525 }
1526
1527 CHKERR printPartitionedProblem(&*out_problem_it, verb);
1528 CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1529
1530 // Inidcate that porble has been build
1533
1535}
#define ProblemManagerFunctionBegin
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ 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.
virtual const Problem_multiIndex * get_problems() const =0
Get the problems object.
virtual MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)=0
clear problem
const double v
phase velocity of light in medium (cm/ns)
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
int & getBuildMoFEM() const
Get flags/semaphores for different stages.
Definition Core.hpp:225
Deprecated interface functions.
MoFEMErrorCode printPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)
MoFEMErrorCode debugPartitionedProblem(const Problem *problem_ptr, int verb=VERBOSE)

◆ buildProblem() [1/2]

MoFEMErrorCode MoFEM::ProblemsManager::buildProblem ( const std::string name,
const bool square_matrix,
int verb = VERBOSE )

build problem data structures

Note
If square_matrix is set to true, that indicate that problem is structurally symmetric, i.e. rows and columns have the same dofs and are indexed in the same way.
Parameters
nameproblem name
square_matrixstructurally symmetric problem
verbverbosity level
Returns
error code
Examples
build_problems.cpp, continuity_check_on_contact_prism_side_ele.cpp, continuity_check_on_skeleton_3d.cpp, forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, hcurl_divergence_operator_2d.cpp, mesh_insert_interface_atom.cpp, nonlinear_dynamics.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, and remove_entities_from_problem_not_partitioned.cpp.

Definition at line 87 of file ProblemsManager.cpp.

89 {
90
91 MoFEM::Interface &m_field = cOre;
93 if (!(cOre.getBuildMoFEM() & (1 << 0)))
94 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
95 if (!(cOre.getBuildMoFEM() & (1 << 1)))
96 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
97 if (!(cOre.getBuildMoFEM() & (1 << 2)))
98 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
99 const Problem *problem_ptr;
100 CHKERR m_field.get_problem(name, &problem_ptr);
101 CHKERR buildProblem(const_cast<Problem *>(problem_ptr), square_matrix, verb);
102 cOre.getBuildMoFEM() |= 1 << 3; // It is assumed that user who uses this
103 // function knows what he is doing
105}
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures

◆ buildProblem() [2/2]

MoFEMErrorCode MoFEM::ProblemsManager::buildProblem ( Problem * problem_ptr,
const bool square_matrix,
int verb = VERBOSE )

build problem data structures

Note
If square_matrix is set to true, that indicate that problem is structurally symmetric, i.e. rows and columns have the same dofs and are indexed in the same way.
Parameters
problempointer
square_matrixstructurally symmetric problem
verbverbosity level
Returns
error code

Definition at line 107 of file ProblemsManager.cpp.

109 {
110 MoFEM::Interface &m_field = cOre;
111 auto dofs_field_ptr = m_field.get_dofs();
112 auto ents_field_ptr = m_field.get_field_ents();
113 auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
115 PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
116
117 // Note: Only allowed changes on problem_ptr structure which not influence
118 // multi-index.
119
120 if (problem_ptr->getBitRefLevel().none()) {
121 SETERRQ1(PETSC_COMM_SELF, 1, "problem <%s> refinement level not set",
122 problem_ptr->getName().c_str());
123 }
124 CHKERR m_field.clear_problem(problem_ptr->getName());
125
126 // zero finite elements
127 problem_ptr->numeredFiniteElementsPtr->clear();
128
129 DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
130 auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
132 for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end(); ++it) {
133
134 const auto uid = (*it)->getLocalUniqueId();
135
136 auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
137 for (auto lo = r.first; lo != r.second; ++lo) {
138
139 if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
140 std::array<bool, 2> row_col = {false, false};
141
142 const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
143 const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
144 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
145
146 // if entity is not problem refinement level
147 if ((fe_bit & prb_mask) != fe_bit)
148 continue;
149 if ((fe_bit & prb_bit).none())
150 continue;
151
152 auto add_to_view = [&](auto &nb_dofs, auto &view, auto rc) {
153 auto dit = nb_dofs->lower_bound(uid);
154 decltype(dit) hi_dit;
155 if (dit != nb_dofs->end()) {
156 hi_dit = nb_dofs->upper_bound(
157 uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
158 view.insert(dit, hi_dit);
159 row_col[rc] = true;
160 }
161 };
162
163 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
164 add_to_view(dofs_field_ptr, dofs_rows, ROW);
165
166 if (!square_matrix)
167 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
168 add_to_view(dofs_field_ptr, dofs_cols, COL);
169
170 if (square_matrix && row_col[ROW])
171 break;
172 else if (row_col[ROW] && row_col[COL])
173 break;
174 }
175 }
176 }
178 };
179
180 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
181
182 // Add row dofs to problem
183 {
184 // zero rows
185 problem_ptr->nbDofsRow = 0;
186 problem_ptr->nbLocDofsRow = 0;
187 problem_ptr->nbGhostDofsRow = 0;
188
189 // add dofs for rows
190 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
191 hi_miit;
192 hi_miit = dofs_rows.get<0>().end();
193
194 int count_dofs = dofs_rows.get<1>().count(true);
195 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
196 boost::shared_ptr<std::vector<NumeredDofEntity>>(
197 new std::vector<NumeredDofEntity>());
198 problem_ptr->getRowDofsSequence()->push_back(dofs_array);
199 dofs_array->reserve(count_dofs);
200 miit = dofs_rows.get<0>().begin();
201 for (; miit != hi_miit; miit++) {
202 if ((*miit)->getActive()) {
203 dofs_array->emplace_back(*miit);
204 dofs_array->back().dofIdx = (problem_ptr->nbDofsRow)++;
205 }
206 }
207 auto hint = problem_ptr->numeredRowDofsPtr->end();
208 for (auto &v : *dofs_array) {
209 hint = problem_ptr->numeredRowDofsPtr->emplace_hint(hint, dofs_array, &v);
210 }
211 }
212
213 // Add col dofs to problem
214 if (!square_matrix) {
215 // zero cols
216 problem_ptr->nbDofsCol = 0;
217 problem_ptr->nbLocDofsCol = 0;
218 problem_ptr->nbGhostDofsCol = 0;
219
220 // add dofs for cols
221 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
222 hi_miit;
223 hi_miit = dofs_cols.get<0>().end();
224
225 int count_dofs = 0;
226 miit = dofs_cols.get<0>().begin();
227 for (; miit != hi_miit; miit++) {
228 if (!(*miit)->getActive()) {
229 continue;
230 }
231 count_dofs++;
232 }
233
234 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
235 boost::shared_ptr<std::vector<NumeredDofEntity>>(
236 new std::vector<NumeredDofEntity>());
237 problem_ptr->getColDofsSequence()->push_back(dofs_array);
238 dofs_array->reserve(count_dofs);
239 miit = dofs_cols.get<0>().begin();
240 for (; miit != hi_miit; miit++) {
241 if (!(*miit)->getActive()) {
242 continue;
243 }
244 dofs_array->emplace_back(*miit);
245 dofs_array->back().dofIdx = problem_ptr->nbDofsCol++;
246 }
247 auto hint = problem_ptr->numeredColDofsPtr->end();
248 for (auto &v : *dofs_array) {
249 hint = problem_ptr->numeredColDofsPtr->emplace_hint(hint, dofs_array, &v);
250 }
251 } else {
252 problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
253 problem_ptr->nbLocDofsCol = problem_ptr->nbLocDofsRow;
254 problem_ptr->nbDofsCol = problem_ptr->nbDofsRow;
255 }
256
257 // job done, some debugging and postprocessing
258 if (verb >= VERBOSE) {
259 MOFEM_LOG("SYNC", Sev::verbose)
260 << problem_ptr->getName() << " Nb. local dofs "
261 << problem_ptr->numeredRowDofsPtr->size() << " by "
262 << problem_ptr->numeredColDofsPtr->size();
263 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
264 }
265
266 if (verb >= NOISY) {
267 MOFEM_LOG("SYNC", Sev::noisy)
268 << "FEs row dofs " << *problem_ptr << " Nb. row dof "
269 << problem_ptr->getNbDofsRow();
270 for (auto &miit : *problem_ptr->numeredRowDofsPtr)
271 MOFEM_LOG("SYNC", Sev::noisy) << *miit;
272
273 MOFEM_LOG("SYNC", Sev::noisy)
274 << "FEs col dofs " << *problem_ptr << " Nb. col dof "
275 << problem_ptr->getNbDofsCol();
276 for (auto &miit : *problem_ptr->numeredColDofsPtr)
277 MOFEM_LOG("SYNC", Sev::noisy) << *miit;
278 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
279 }
280
281 cOre.getBuildMoFEM() |= Core::BUILD_PROBLEM; // It is assumed that user who
282 // uses this function knows
283 // what he is doing
284
285 PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
286
288}
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
@ NOISY
@ VERBOSE
@ COL
@ ROW
#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 MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getLocalUniqueId > >, ordered_non_unique< const_mem_fun< DofEntity, char, &DofEntity::getActive > > > > DofEntity_multiIndex_active_view
multi-index view on DofEntity activity
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
virtual MoFEMErrorCode get_ents_elements_adjacency(const FieldEntityEntFiniteElementAdjacencyMap_multiIndex **dofs_elements_adjacency) const =0
Get the dofs elements adjacency object.
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
#define MOFEM_LOG(channel, severity)
Log.
uint128_t UId
Unique Id.
Definition Types.hpp:31
int r
Definition sdf.py:8
PetscLogEvent MOFEM_EVENT_ProblemsManager

◆ buildProblemOnDistributedMesh() [1/2]

MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh ( const std::string name,
const bool square_matrix,
int verb = VERBOSE )

build problem data structures, assuming that mesh is distributed (collective)

Mesh is distributed, that means that each processor keeps only own part of the mesh and shared entities.

Collective - need to be run on all processors in communicator, i.e. each function has to call this function.

Examples
nonlinear_dynamics.cpp, and remove_entities_from_problem.cpp.

Definition at line 290 of file ProblemsManager.cpp.

291 {
292 MoFEM::Interface &m_field = cOre;
294
296 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
297 if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
298 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
300 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
301
302 const Problem *problem_ptr;
303 CHKERR m_field.get_problem(name, &problem_ptr);
304 CHKERR buildProblemOnDistributedMesh(const_cast<Problem *>(problem_ptr),
305 square_matrix, verb);
306
309
311}
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)

◆ buildProblemOnDistributedMesh() [2/2]

MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh ( Problem * problem_ptr,
const bool square_matrix = true,
int verb = VERBOSE )

build problem data structures, assuming that mesh is distributed (collective)

Mesh is distributed, that means that each processor keeps only own part of the mesh and shared entities.

Collective - need to be run on all processors in communicator, i.e. each function has to call this function.

Definition at line 313 of file ProblemsManager.cpp.

314 {
315 MoFEM::Interface &m_field = cOre;
316 auto fields_ptr = m_field.get_fields();
317 auto fe_ptr = m_field.get_finite_elements();
318 auto dofs_field_ptr = m_field.get_dofs();
320 PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
321
322 // clear data structures
323 CHKERR m_field.clear_problem(problem_ptr->getName());
324
326
327 if (problem_ptr->getBitRefLevel().none())
328 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
329 "problem <%s> refinement level not set",
330 problem_ptr->getName().c_str());
331
332 int loop_size = 2;
333 if (square_matrix) {
334 loop_size = 1;
335 problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
336 } else if (problem_ptr->numeredColDofsPtr == problem_ptr->numeredRowDofsPtr) {
337 problem_ptr->numeredColDofsPtr =
338 boost::shared_ptr<NumeredDofEntity_multiIndex>(
340 }
341
342 const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
343 const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
344
345 // // get rows and cols dofs view based on data on elements
346 DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
347
348 // Add DOFs to problem by visiting all elements and adding DOFs from
349 // elements to the problem
350 if (buildProblemFromFields == PETSC_FALSE) {
351
352 auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
353 auto ents_field_ptr = m_field.get_field_ents();
354 auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
356 for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end();
357 ++it) {
358
359 const auto uid = (*it)->getLocalUniqueId();
360
361 auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
362 for (auto lo = r.first; lo != r.second; ++lo) {
363
364 if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
365 std::array<bool, 2> row_col = {false, false};
366
367 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
368
369 // if entity is not problem refinement level
370 if ((fe_bit & prb_bit).any() && (fe_bit & prb_mask) == fe_bit) {
371
372 auto add_to_view = [&](auto dofs, auto &view, auto rc) {
373 auto dit = dofs->lower_bound(uid);
374 decltype(dit) hi_dit;
375 if (dit != dofs->end()) {
376 hi_dit = dofs->upper_bound(
377 uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
378 view.insert(dit, hi_dit);
379 row_col[rc] = true;
380 }
381 };
382
383 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
384 add_to_view(dofs_field_ptr, dofs_rows, ROW);
385
386 if (!square_matrix)
387 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
388 add_to_view(dofs_field_ptr, dofs_cols, COL);
389
390 if (square_matrix && row_col[ROW])
391 break;
392 else if (row_col[ROW] && row_col[COL])
393 break;
394 }
395 }
396 }
397 }
399 };
400
401 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
402 }
403
404 // Add DOFS to the problem by searching all the filedes, and adding to
405 // problem owned or shared DOFs
406 if (buildProblemFromFields == PETSC_TRUE) {
407 // Get fields IDs on elements
408 BitFieldId fields_ids_row, fields_ids_col;
409 for (auto fit = fe_ptr->begin(); fit != fe_ptr->end(); fit++) {
410 if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
411 fields_ids_row |= fit->get()->getBitFieldIdRow();
412 fields_ids_col |= fit->get()->getBitFieldIdCol();
413 }
414 }
415 // Get fields DOFs
416 for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
417 if ((fit->get()->getId() & (fields_ids_row | fields_ids_col)).any()) {
418
419 auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(
420 FieldEntity::getLoBitNumberUId((*fit)->getBitNumber()));
421 auto hi_dit = dofs_field_ptr->get<Unique_mi_tag>().upper_bound(
422 FieldEntity::getHiBitNumberUId((*fit)->getBitNumber()));
423
424 for (; dit != hi_dit; dit++) {
425
426 const int owner_proc = dit->get()->getOwnerProc();
427 if (owner_proc != m_field.get_comm_rank()) {
428 const unsigned char pstatus = dit->get()->getPStatus();
429 if (pstatus == 0) {
430 continue;
431 }
432 }
433
434 const auto &dof_bit = (*dit)->getBitRefLevel();
435 // if entity is not problem refinement level
436 if ((dof_bit & prb_bit).any() && (dof_bit & prb_mask) == dof_bit) {
437
438 if ((fit->get()->getId() & fields_ids_row).any()) {
439 dofs_rows.insert(*dit);
440 }
441 if (!square_matrix) {
442 if ((fit->get()->getId() & fields_ids_col).any()) {
443 dofs_cols.insert(*dit);
444 }
445 }
446 }
447 }
448 }
449 }
450 }
451
453 // Get fields IDs on elements
454 BitFieldId fields_ids_row, fields_ids_col;
455 BitFieldId *fields_ids[2] = {&fields_ids_row, &fields_ids_col};
456 for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
457 fit != fe_ptr->end(); fit++) {
458 if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
459 fields_ids_row |= fit->get()->getBitFieldIdRow();
460 fields_ids_col |= fit->get()->getBitFieldIdCol();
461 }
462 }
463
464 DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
465 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ++ss) {
466 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
467 hi_miit;
468 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
469 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
470 Range ents_to_synchronise;
471 for (; miit != hi_miit; ++miit) {
472 if (miit->get()->getEntDofIdx() != 0)
473 continue;
474 ents_to_synchronise.insert(miit->get()->getEnt());
475 }
476 Range tmp_ents = ents_to_synchronise;
477 CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
478 ents_to_synchronise, nullptr, verb);
479 ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
480 for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
481 if ((fit->get()->getId() & *fields_ids[ss]).any()) {
482 const auto bit_number = (*fit)->getBitNumber();
483 for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
484 pit != ents_to_synchronise.pair_end(); ++pit) {
485 const auto f = pit->first;
486 const auto s = pit->second;
487 const auto lo_uid =
489 const auto hi_uid =
491
492 auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
493 auto hi_dit =
494 dofs_field_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
495
496 dofs_ptr[ss]->insert(dit, hi_dit);
497 }
498 }
499 }
500 }
501 }
502
503 // add dofs for rows and cols and set ownership
504 DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
505 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
506 problem_ptr->numeredRowDofsPtr, problem_ptr->numeredColDofsPtr};
507 int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
508 int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
509 &problem_ptr->nbLocDofsCol};
510 int *ghost_nbdof_ptr[] = {&problem_ptr->nbGhostDofsRow,
511 &problem_ptr->nbGhostDofsCol};
512 for (int ss = 0; ss < 2; ss++) {
513 *(nbdof_ptr[ss]) = 0;
514 *(local_nbdof_ptr[ss]) = 0;
515 *(ghost_nbdof_ptr[ss]) = 0;
516 }
517
518 // Loop over dofs on rows and columns and add to multi-indices in dofs
519 // problem structure, set partition for each dof
520 int nb_local_dofs[] = {0, 0};
521 for (int ss = 0; ss < loop_size; ss++) {
522 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
523 hi_miit;
524 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
525 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
526 for (; miit != hi_miit; miit++) {
527 int owner_proc = (*miit)->getOwnerProc();
528 if (owner_proc == m_field.get_comm_rank()) {
529 nb_local_dofs[ss]++;
530 }
531 }
532 }
533
534 // get layout
535 int start_ranges[2], end_ranges[2];
536 for (int ss = 0; ss != loop_size; ss++) {
537 PetscLayout layout;
538 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
539 CHKERR PetscLayoutSetBlockSize(layout, 1);
540 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
541 CHKERR PetscLayoutSetUp(layout);
542 CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
543 CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
544 CHKERR PetscLayoutDestroy(&layout);
545 }
546 if (square_matrix) {
547 nbdof_ptr[1] = nbdof_ptr[0];
548 nb_local_dofs[1] = nb_local_dofs[0];
549 start_ranges[1] = start_ranges[0];
550 end_ranges[1] = end_ranges[0];
551 }
552
553 // if(sizeof(UId) != SIZEOFUID) {
554 // SETERRQ2(
555 // PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,
556 // "check size of UId, size of UId is %u != %u",
557 // sizeof(UId),SIZEOFUID
558 // );
559 // }
560
561 // set local and global indices on own dofs
562 const size_t idx_data_type_size = sizeof(IdxDataType);
563 const size_t data_block_size = idx_data_type_size / sizeof(int);
564
565 if (sizeof(IdxDataType) % sizeof(int)) {
566 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
567 }
568 std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
569 m_field.get_comm_size()),
570 ids_data_packed_cols(m_field.get_comm_size());
571
572 // Loop over dofs on this processor and prepare those dofs to send on
573 // another proc
574 for (int ss = 0; ss != loop_size; ++ss) {
575
576 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
577 hi_miit;
578 hi_miit = dofs_ptr[ss]->get<0>().end();
579
580 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
581 boost::shared_ptr<std::vector<NumeredDofEntity>>(
582 new std::vector<NumeredDofEntity>());
583 int nb_dofs_to_add = 0;
584 miit = dofs_ptr[ss]->get<0>().begin();
585 for (; miit != hi_miit; ++miit) {
586 // Only set global idx for dofs on this processor part
587 if (!(miit->get()->getActive()))
588 continue;
589 ++nb_dofs_to_add;
590 }
591 dofs_array->reserve(nb_dofs_to_add);
592 if (ss == 0) {
593 problem_ptr->getRowDofsSequence()->push_back(dofs_array);
594 } else {
595 problem_ptr->getColDofsSequence()->push_back(dofs_array);
596 }
597
598 int &local_idx = *local_nbdof_ptr[ss];
599 miit = dofs_ptr[ss]->get<0>().begin();
600 for (; miit != hi_miit; ++miit) {
601
602 // Only set global idx for dofs on this processor part
603 if (!(miit->get()->getActive()))
604 continue;
605
606 dofs_array->emplace_back(*miit);
607
608 int owner_proc = dofs_array->back().getOwnerProc();
609 if (owner_proc < 0) {
610 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
611 "data inconsistency");
612 }
613
614 if (owner_proc != m_field.get_comm_rank()) {
615 dofs_array->back().pArt = owner_proc;
616 dofs_array->back().dofIdx = -1;
617 dofs_array->back().petscGloablDofIdx = -1;
618 dofs_array->back().petscLocalDofIdx = -1;
619 } else {
620
621 // Set part and indexes
622 int glob_idx = start_ranges[ss] + local_idx;
623 dofs_array->back().pArt = owner_proc;
624 dofs_array->back().dofIdx = glob_idx;
625 dofs_array->back().petscGloablDofIdx = glob_idx;
626 dofs_array->back().petscLocalDofIdx = local_idx;
627 local_idx++;
628
629 unsigned char pstatus = dofs_array->back().getPStatus();
630 // check id dof is shared, if that is a case global idx added to data
631 // structure and is sended to other processors
632 if (pstatus > 0) {
633
634 for (int proc = 0;
635 proc < MAX_SHARING_PROCS &&
636 -1 != dofs_array->back().getSharingProcsPtr()[proc];
637 proc++) {
638 // make it different for rows and columns when needed
639 if (ss == 0) {
640 ids_data_packed_rows[dofs_array->back()
641 .getSharingProcsPtr()[proc]]
642 .emplace_back(dofs_array->back().getGlobalUniqueId(),
643 glob_idx);
644 } else {
645 ids_data_packed_cols[dofs_array->back()
646 .getSharingProcsPtr()[proc]]
647 .emplace_back(dofs_array->back().getGlobalUniqueId(),
648 glob_idx);
649 }
650 if (!(pstatus & PSTATUS_MULTISHARED)) {
651 break;
652 }
653 }
654 }
655 }
656 }
657
658 auto hint = numered_dofs_ptr[ss]->end();
659 for (auto &v : *dofs_array)
660 hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &v);
661 }
662 if (square_matrix) {
663 local_nbdof_ptr[1] = local_nbdof_ptr[0];
664 }
665
666 int nsends_rows = 0, nsends_cols = 0;
667 // Non zero lengths[i] represent a message to i of length lengths[i].
668 std::vector<int> lengths_rows(m_field.get_comm_size()),
669 lengths_cols(m_field.get_comm_size());
670 lengths_rows.clear();
671 lengths_cols.clear();
672 for (int proc = 0; proc < m_field.get_comm_size(); proc++) {
673 lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
674 lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
675 if (!ids_data_packed_rows[proc].empty())
676 nsends_rows++;
677 if (!ids_data_packed_cols[proc].empty())
678 nsends_cols++;
679 }
680
681 MPI_Status *status;
682 CHKERR PetscMalloc1(m_field.get_comm_size(), &status);
683
684 // Do rows
685 int nrecvs_rows; // number of messages received
686 int *onodes_rows; // list of node-ids from which messages are expected
687 int *olengths_rows; // corresponding message lengths
688 int **rbuf_row; // must bee freed by user
689
690 // make sure it is a PETSc comm
691 MPI_Comm comm;
692 CHKERR PetscCommDuplicate(m_field.get_comm(), &comm, NULL);
693
694 {
695
696 // rows
697
698 // Computes the number of messages a node expects to receive
699 CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_rows[0],
700 &nrecvs_rows);
701 // std::cerr << nrecvs_rows << std::endl;
702
703 // Computes info about messages that a MPI-node will receive, including
704 // (from-id,length) pairs for each message.
705 CHKERR PetscGatherMessageLengths(comm, nsends_rows, nrecvs_rows,
706 &lengths_rows[0], &onodes_rows,
707 &olengths_rows);
708
709 // Gets a unique new tag from a PETSc communicator. All processors that
710 // share the communicator MUST call this routine EXACTLY the same number
711 // of times. This tag should only be used with the current objects
712 // communicator; do NOT use it with any other MPI communicator.
713 int tag_row;
714 CHKERR PetscCommGetNewTag(comm, &tag_row);
715
716 // Allocate a buffer sufficient to hold messages of size specified in
717 // olengths. And post Irecvs on these buffers using node info from onodes
718 MPI_Request *r_waits_row; // must bee freed by user
719 // rbuf has a pointers to messeges. It has size of of nrecvs (number of
720 // messages) +1. In the first index a block is allocated,
721 // such that rbuf[i] = rbuf[i-1]+olengths[i-1].
722
723 CHKERR PetscPostIrecvInt(comm, tag_row, nrecvs_rows, onodes_rows,
724 olengths_rows, &rbuf_row, &r_waits_row);
725 CHKERR PetscFree(onodes_rows);
726
727 MPI_Request *s_waits_row; // status of sens messages
728 CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
729
730 // Send messeges
731 for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
732 if (!lengths_rows[proc])
733 continue; // no message to send to this proc
734 CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0], // buffer to send
735 lengths_rows[proc], // message length
736 MPIU_INT, proc, // to proc
737 tag_row, comm, s_waits_row + kk);
738 kk++;
739 }
740
741 if (nrecvs_rows) {
742 CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
743 }
744 if (nsends_rows) {
745 CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
746 }
747
748 CHKERR PetscFree(r_waits_row);
749 CHKERR PetscFree(s_waits_row);
750 }
751
752 // cols
753 int nrecvs_cols = nrecvs_rows;
754 int *olengths_cols = olengths_rows;
755 PetscInt **rbuf_col = rbuf_row;
756 if (!square_matrix) {
757
758 // Computes the number of messages a node expects to receive
759 CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_cols[0],
760 &nrecvs_cols);
761
762 // Computes info about messages that a MPI-node will receive, including
763 // (from-id,length) pairs for each message.
764 int *onodes_cols;
765 CHKERR PetscGatherMessageLengths(comm, nsends_cols, nrecvs_cols,
766 &lengths_cols[0], &onodes_cols,
767 &olengths_cols);
768
769 // Gets a unique new tag from a PETSc communicator.
770 int tag_col;
771 CHKERR PetscCommGetNewTag(comm, &tag_col);
772
773 // Allocate a buffer sufficient to hold messages of size specified in
774 // olengths. And post Irecvs on these buffers using node info from onodes
775 MPI_Request *r_waits_col; // must bee freed by user
776 CHKERR PetscPostIrecvInt(comm, tag_col, nrecvs_cols, onodes_cols,
777 olengths_cols, &rbuf_col, &r_waits_col);
778 CHKERR PetscFree(onodes_cols);
779
780 MPI_Request *s_waits_col; // status of sens messages
781 CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
782
783 // Send messages
784 for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
785 if (!lengths_cols[proc])
786 continue; // no message to send to this proc
787 CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0], // buffer to send
788 lengths_cols[proc], // message length
789 MPIU_INT, proc, // to proc
790 tag_col, comm, s_waits_col + kk);
791 kk++;
792 }
793
794 if (nrecvs_cols) {
795 CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
796 }
797 if (nsends_cols) {
798 CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
799 }
800
801 CHKERR PetscFree(r_waits_col);
802 CHKERR PetscFree(s_waits_col);
803 }
804
805 CHKERR PetscCommDestroy(&comm);
806 CHKERR PetscFree(status);
807
808 DofEntity_multiIndex_global_uid_view dofs_glob_uid_view;
809 auto hint = dofs_glob_uid_view.begin();
810 for (auto dof : *m_field.get_dofs())
811 dofs_glob_uid_view.emplace_hint(hint, dof);
812
813 // set values received from other processors
814 for (int ss = 0; ss != loop_size; ++ss) {
815
816 int nrecvs;
817 int *olengths;
818 int **data_procs;
819 if (ss == 0) {
820 nrecvs = nrecvs_rows;
821 olengths = olengths_rows;
822 data_procs = rbuf_row;
823 } else {
824 nrecvs = nrecvs_cols;
825 olengths = olengths_cols;
826 data_procs = rbuf_col;
827 }
828
829 UId uid;
830 for (int kk = 0; kk != nrecvs; ++kk) {
831 int len = olengths[kk];
832 int *data_from_proc = data_procs[kk];
833 for (int dd = 0; dd < len; dd += data_block_size) {
834 uid = IdxDataTypePtr(&data_from_proc[dd]).getUId();
835 auto ddit = dofs_glob_uid_view.find(uid);
836
837 if (PetscUnlikely(ddit == dofs_glob_uid_view.end())) {
838
839#ifndef NDEBUG
840 auto ents_field_ptr = m_field.get_field_ents();
841 MOFEM_LOG("SELF", Sev::error)
842 << "DOF is shared or multishared between processors. For example "
843 "if order of field on given entity is set in inconsistently, "
844 "has different value on two processor, error such as this is "
845 "triggered";
846
847 MOFEM_LOG("SELF", Sev::error) << "UId " << uid << " is not found";
848 const auto owner_proc = FieldEntity::getOwnerFromUniqueId(uid);
849 MOFEM_LOG("SELF", Sev::error)
850 << "Problematic UId owner proc is " << owner_proc;
851 const auto uid_handle = FieldEntity::getHandleFromUniqueId(uid);
852 MOFEM_LOG("SELF", Sev::error)
853 << "Problematic UId entity owning processor handle is "
854 << uid_handle << " entity type "
855 << moab::CN::EntityTypeName(type_from_handle(uid_handle));
856 const auto uid_bit_number =
858 MOFEM_LOG("SELF", Sev::error)
859 << "Problematic UId field is "
860 << m_field.get_field_name(
861 field_bit_from_bit_number(uid_bit_number));
862
864 field_view.insert(ents_field_ptr->begin(), ents_field_ptr->end());
865 auto fe_it = field_view.find(FieldEntity::getGlobalUniqueIdCalculate(
866 owner_proc, uid_bit_number, uid_handle));
867 if (fe_it == field_view.end()) {
868 MOFEM_LOG("SELF", Sev::error)
869 << "Also, no field entity in database for given global UId";
870 } else {
871 MOFEM_LOG("SELF", Sev::error) << "Field entity in databse exist "
872 "(but have no DOF wih give UId";
873 MOFEM_LOG("SELF", Sev::error) << **fe_it;
874
875 // Save file with missing entity
876 auto error_file_name =
877 "error_with_missing_entity_" +
878 boost::lexical_cast<std::string>(m_field.get_comm_rank()) +
879 ".vtk";
880 MOFEM_LOG("SELF", Sev::error)
881 << "Look to file < " << error_file_name
882 << " > it contains entity with missing DOF.";
883
884 auto tmp_msh = get_temp_meshset_ptr(m_field.get_moab());
885 const auto local_fe_ent = (*fe_it)->getEnt();
886 CHKERR m_field.get_moab().add_entities(*tmp_msh, &local_fe_ent, 1);
887 CHKERR m_field.get_moab().write_file(error_file_name.c_str(), "VTK",
888 "", tmp_msh->get_ptr(), 1);
889 }
890#endif
891
892 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
893 "DOF with global UId not found (Compile code in Debug to "
894 "learn more about problem");
895 }
896
897 auto dit = numered_dofs_ptr[ss]->find((*ddit)->getLocalUniqueId());
898
899 if (dit != numered_dofs_ptr[ss]->end()) {
900
901 int global_idx = IdxDataTypePtr(&data_from_proc[dd]).getDofIdx();
902#ifndef NDEBUG
903 if (PetscUnlikely(global_idx < 0))
904 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
905 "received negative dof");
906#endif
907 bool success;
908 success = numered_dofs_ptr[ss]->modify(
909 dit, NumeredDofEntity_mofem_index_change(global_idx));
910
911#ifndef NDEBUG
912 if (PetscUnlikely(!success))
913 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
914 "modification unsuccessful");
915#endif
916 success = numered_dofs_ptr[ss]->modify(
917 dit, NumeredDofEntity_part_and_glob_idx_change((*dit)->getPart(),
918 global_idx));
919#ifndef NDEBUG
920 if (PetscUnlikely(!success))
921 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
922 "modification unsuccessful");
923#endif
924 };
925
926#ifndef NDEBUG
927 if (PetscUnlikely(ddit->get()->getPStatus() == 0)) {
928
929 // Dof is shared on this processor, however there is no element
930 // which have this dof. If DOF is not shared and received from other
931 // processor, but not marked as a shared on other that means that is
932 // data inconstancy and error should be thrown.
933
934 std::ostringstream zz;
935 zz << **ddit << std::endl;
936 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
937 "data inconsistency, dofs is not shared, but received "
938 "from other proc\n"
939 "%s",
940 zz.str().c_str());
941 }
942#endif
943 }
944 }
945 }
946
947 if (square_matrix) {
948 (problem_ptr->nbDofsCol) = (problem_ptr->nbDofsRow);
949 (problem_ptr->nbLocDofsCol) = (problem_ptr->nbLocDofsRow);
950 }
951
952 CHKERR PetscFree(olengths_rows);
953 CHKERR PetscFree(rbuf_row[0]);
954 CHKERR PetscFree(rbuf_row);
955 if (!square_matrix) {
956 CHKERR PetscFree(olengths_cols);
957 CHKERR PetscFree(rbuf_col[0]);
958 CHKERR PetscFree(rbuf_col);
959 }
960
961 if (square_matrix) {
962 if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
963 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
964 "data inconsistency for square_matrix %d!=%d",
965 numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
966 }
967 if (problem_ptr->numeredRowDofsPtr != problem_ptr->numeredColDofsPtr) {
968 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
969 "data inconsistency for square_matrix");
970 }
971 }
972
973 CHKERR printPartitionedProblem(problem_ptr, verb);
974 CHKERR debugPartitionedProblem(problem_ptr, verb);
975
976 PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
977
979}
@ MOFEM_NOT_FOUND
Definition definitions.h:33
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.
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< const_mem_fun< DofEntity, UId, &DofEntity::getGlobalUniqueId > > > > DofEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< FieldEntity, UId, &FieldEntity::getGlobalUniqueId > > > > FieldEntity_multiIndex_global_uid_view
multi-index view on DofEntity by uid
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition Types.hpp:42
auto type_from_handle(const EntityHandle h)
get type from entity handle
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual std::string get_field_name(const BitFieldId id) const =0
get field name from id
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 auto getOwnerFromUniqueId(const UId uid)
UId getGlobalUniqueIdCalculate() const
Calculate global UId.
PetscBool synchroniseProblemEntities
DOFs in fields, not from DOFs on elements.
MoFEMErrorCode getOptions()
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

◆ buildSubProblem()

MoFEMErrorCode MoFEM::ProblemsManager::buildSubProblem ( const std::string out_name,
const std::vector< std::string > & fields_row,
const std::vector< std::string > & fields_col,
const std::string main_problem,
const bool square_matrix = true,
const map< std::string, boost::shared_ptr< Range > > * entityMapRow = nullptr,
const map< std::string, boost::shared_ptr< Range > > * entityMapCol = nullptr,
int verb = VERBOSE )

build sub problem

Parameters
out_nameproblem
fields_rowvector of fields composing problem
fields_colvector of fields composing problem
main_problemmain problem
Returns
error code
Examples
free_surface.cpp, and remove_entities_from_problem.cpp.

Definition at line 981 of file ProblemsManager.cpp.

992 {
993 MoFEM::Interface &m_field = cOre;
994 auto problems_ptr = m_field.get_problems();
996
997 CHKERR m_field.clear_problem(out_name);
998
999 // get reference to all problems
1000 using ProblemByName = decltype(problems_ptr->get<Problem_mi_tag>());
1001 auto &problems_by_name =
1002 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1003
1004 // get iterators to out problem, i.e. build problem
1005 auto out_problem_it = problems_by_name.find(out_name);
1006 if (out_problem_it == problems_by_name.end()) {
1007 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1008 "subproblem with name < %s > not defined (top tip check spelling)",
1009 out_name.c_str());
1010 }
1011 // get iterator to main problem, i.e. out problem is subproblem of main
1012 // problem
1013 auto main_problem_it = problems_by_name.find(main_problem);
1014 if (main_problem_it == problems_by_name.end()) {
1015 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1016 "problem of subproblem with name < %s > not defined (top tip "
1017 "check spelling)",
1018 main_problem.c_str());
1019 }
1020
1021 // get dofs for row & columns for out problem,
1022 boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1023 out_problem_it->numeredRowDofsPtr, out_problem_it->numeredColDofsPtr};
1024 // get dofs for row & columns for main problem
1025 boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1026 main_problem_it->numeredRowDofsPtr, main_problem_it->numeredColDofsPtr};
1027 // get local indices counter
1028 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1029 &out_problem_it->nbLocDofsCol};
1030 // get global indices counter
1031 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1032
1033 // set number of ghost nodes to zero
1034 {
1035 out_problem_it->nbGhostDofsRow = 0;
1036 out_problem_it->nbGhostDofsCol = 0;
1037 }
1038
1039 // put rows & columns field names in array
1040 std::vector<std::string> fields[] = {fields_row, fields_col};
1041 const map<std::string, boost::shared_ptr<Range>> *entityMap[] = {
1042 entityMapRow, entityMapCol};
1043
1044 // make data structure fos sub-problem data
1045 out_problem_it->subProblemData =
1046 boost::make_shared<Problem::SubProblemData>();
1047
1048 // Loop over rows and columns
1049 for (int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1050
1051 // reset dofs and columns counters
1052 (*nb_local_dofs[ss]) = 0;
1053 (*nb_dofs[ss]) = 0;
1054 // clear arrays
1055 out_problem_dofs[ss]->clear();
1056
1057 // If DOFs are cleared clear finite elements too.
1058 out_problem_it->numeredFiniteElementsPtr->clear();
1059
1060 // get dofs by field name and insert them in out problem multi-indices
1061 for (auto field : fields[ss]) {
1062
1063 // Following reserve memory in sequences, only two allocations are here,
1064 // once for array of objects, next for array of shared pointers
1065
1066 // aliased sequence of pointer is killed with element
1067 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1068 boost::make_shared<std::vector<NumeredDofEntity>>();
1069 // reserve memory for field dofs
1070 if (!ss)
1071 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1072 else
1073 out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1074
1075 // create elements objects
1076 auto bit_number = m_field.get_field_bit_number(field);
1077
1078 auto add_dit_to_dofs_array = [&](auto &dit) {
1079 if (dit->get()->getPetscGlobalDofIdx() >= 0)
1080 dofs_array->emplace_back(
1081 dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1082 dit->get()->getPetscGlobalDofIdx(),
1083 dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1084 };
1085
1086 auto get_dafult_dof_range = [&]() {
1087 auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1088 FieldEntity::getLoBitNumberUId(bit_number));
1089 auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1090 FieldEntity::getHiBitNumberUId(bit_number));
1091 return std::make_pair(dit, hi_dit);
1092 };
1093
1094 auto check = [&](auto &field) -> boost::shared_ptr<Range> {
1095 if (entityMap[ss]) {
1096 auto mit = entityMap[ss]->find(field);
1097 if (mit != entityMap[ss]->end()) {
1098 return mit->second;
1099 } else {
1100 return nullptr;
1101 }
1102 } else {
1103 return nullptr;
1104 }
1105 };
1106
1107 if (auto r_ptr = check(field)) {
1108 for (auto p = r_ptr->pair_begin(); p != r_ptr->pair_end(); ++p) {
1109 const auto lo_ent = p->first;
1110 const auto hi_ent = p->second;
1111 auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1112 DofEntity::getLoFieldEntityUId(bit_number, lo_ent));
1113 auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1114 DofEntity::getHiFieldEntityUId(bit_number, hi_ent));
1115 dofs_array->reserve(std::distance(dit, hi_dit));
1116 for (; dit != hi_dit; dit++) {
1117 add_dit_to_dofs_array(dit);
1118 }
1119 }
1120 } else {
1121 auto [dit, hi_dit] = get_dafult_dof_range();
1122 dofs_array->reserve(std::distance(dit, hi_dit));
1123 for (; dit != hi_dit; dit++)
1124 add_dit_to_dofs_array(dit);
1125 }
1126
1127 // fill multi-index
1128 auto hint = out_problem_dofs[ss]->end();
1129 for (auto &v : *dofs_array)
1130 hint = out_problem_dofs[ss]->emplace_hint(hint, dofs_array, &v);
1131 }
1132 // Set local indexes
1133 {
1134 auto dit = out_problem_dofs[ss]->get<Idx_mi_tag>().begin();
1135 auto hi_dit = out_problem_dofs[ss]->get<Idx_mi_tag>().end();
1136 for (; dit != hi_dit; dit++) {
1137 int idx = -1; // if dof is not part of partition, set local index to -1
1138 if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1139 idx = (*nb_local_dofs[ss])++;
1140 }
1141 bool success = out_problem_dofs[ss]->modify(
1142 out_problem_dofs[ss]->project<0>(dit),
1143 NumeredDofEntity_local_idx_change(idx));
1144 if (!success) {
1145 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1146 "operation unsuccessful");
1147 }
1148 };
1149 }
1150 // Set global indexes, compress global indices
1151 {
1152 auto dit =
1153 out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1154 auto hi_dit =
1155 out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(
1156 out_problem_dofs[ss]->size());
1157 const int nb = std::distance(dit, hi_dit);
1158 // get main problem global indices
1159 std::vector<int> main_indices(nb);
1160 for (auto it = main_indices.begin(); dit != hi_dit; dit++, it++) {
1161 *it = dit->get()->getPetscGlobalDofIdx();
1162 }
1163 // create is with global dofs
1164 IS is;
1165 CHKERR ISCreateGeneral(m_field.get_comm(), nb, &*main_indices.begin(),
1166 PETSC_USE_POINTER, &is);
1167 // create map form main problem global indices to out problem global
1168 // indices
1169 AO ao;
1170 CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1171 if (ss == 0) {
1172 IS is_dup;
1173 CHKERR ISDuplicate(is, &is_dup);
1174 out_problem_it->getSubData()->rowIs = SmartPetscObj<IS>(is_dup, false);
1175 out_problem_it->getSubData()->rowMap = SmartPetscObj<AO>(ao, true);
1176 } else {
1177 IS is_dup;
1178 CHKERR ISDuplicate(is, &is_dup);
1179 out_problem_it->getSubData()->colIs = SmartPetscObj<IS>(is_dup, false);
1180 out_problem_it->getSubData()->colMap = SmartPetscObj<AO>(ao, true);
1181 }
1182 CHKERR AOApplicationToPetscIS(ao, is);
1183 // set global number of DOFs
1184 CHKERR ISGetSize(is, nb_dofs[ss]);
1185 CHKERR ISDestroy(&is);
1186 // set out problem global indices after applying map
1187 dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(0);
1188 for (std::vector<int>::iterator it = main_indices.begin(); dit != hi_dit;
1189 dit++, it++) {
1190 bool success = out_problem_dofs[ss]->modify(
1191 out_problem_dofs[ss]->project<0>(dit),
1192 NumeredDofEntity_part_and_all_indices_change(
1193 dit->get()->getPart(), *it, *it,
1194 dit->get()->getPetscLocalDofIdx()));
1195 if (!success) {
1196 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1197 "operation unsuccessful");
1198 }
1199 }
1200 // set global indices to nodes not on this part
1201 {
1202 NumeredDofEntityByLocalIdx::iterator dit =
1203 out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1204 NumeredDofEntityByLocalIdx::iterator hi_dit =
1205 out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().upper_bound(-1);
1206 const int nb = std::distance(dit, hi_dit);
1207 std::vector<int> main_indices_non_local(nb);
1208 for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1209 dit++, it++) {
1210 *it = dit->get()->getPetscGlobalDofIdx();
1211 }
1212 IS is;
1213 CHKERR ISCreateGeneral(m_field.get_comm(), nb,
1214 &*main_indices_non_local.begin(),
1215 PETSC_USE_POINTER, &is);
1216 CHKERR AOApplicationToPetscIS(ao, is);
1217 CHKERR ISDestroy(&is);
1218 dit = out_problem_dofs[ss]->get<PetscLocalIdx_mi_tag>().lower_bound(-1);
1219 for (auto it = main_indices_non_local.begin(); dit != hi_dit;
1220 dit++, it++) {
1221 bool success = out_problem_dofs[ss]->modify(
1222 out_problem_dofs[ss]->project<0>(dit),
1223 NumeredDofEntity_part_and_all_indices_change(
1224 dit->get()->getPart(), dit->get()->getDofIdx(), *it,
1225 dit->get()->getPetscLocalDofIdx()));
1226 if (!success) {
1227 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1228 "operation unsuccessful");
1229 }
1230 }
1231 }
1232 CHKERR AODestroy(&ao);
1233 }
1234 }
1235
1236 if (square_matrix) {
1237 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1238 out_problem_it->nbLocDofsCol = out_problem_it->nbLocDofsRow;
1239 out_problem_it->nbDofsCol = out_problem_it->nbDofsRow;
1240 out_problem_it->getSubData()->colIs = out_problem_it->getSubData()->rowIs;
1241 out_problem_it->getSubData()->colMap = out_problem_it->getSubData()->rowMap;
1242 }
1243
1244 CHKERR printPartitionedProblem(&*out_problem_it, verb);
1245 CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1246
1248}
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)

◆ getFEMeshset()

MoFEMErrorCode MoFEM::ProblemsManager::getFEMeshset ( const std::string prb_name,
const std::string & fe_name,
EntityHandle * meshset ) const

create add entities of finite element in the problem

\create add entities of finite element in the problem

Note
Meshset entity has to be crated
Meshset entity has to be crated
Parameters
prb_namename of the problem
fe_namename of the entity
meshsetpointer meshset handle
Returns
MoFEMErrorCode

Definition at line 2559 of file ProblemsManager.cpp.

2561 {
2562 MoFEM::Interface &m_field = cOre;
2563 const Problem *problem_ptr;
2564 const FiniteElement_multiIndex *fes_ptr;
2566
2567 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2568 CHKERR m_field.get_problem(prb_name, &problem_ptr);
2569 CHKERR m_field.get_finite_elements(&fes_ptr);
2570
2571 auto fe_miit = fes_ptr->get<FiniteElement_name_mi_tag>().find(fe_name);
2572 if (fe_miit != fes_ptr->get<FiniteElement_name_mi_tag>().end()) {
2573 auto fit =
2574 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
2576 0, (*fe_miit)->getFEUId()));
2577 auto hi_fe_it =
2578 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
2580 get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
2581 std::vector<EntityHandle> fe_vec;
2582 fe_vec.reserve(std::distance(fit, hi_fe_it));
2583 for (; fit != hi_fe_it; fit++)
2584 fe_vec.push_back(fit->get()->getEnt());
2585 CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2586 fe_vec.size());
2587 }
2588
2590}
multi_index_container< boost::shared_ptr< FiniteElement >, indexed_by< hashed_unique< tag< FiniteElement_Meshset_mi_tag >, member< FiniteElement, EntityHandle, &FiniteElement::meshset > >, hashed_unique< tag< BitFEId_mi_tag >, const_mem_fun< FiniteElement, BitFEId, &FiniteElement::getId >, HashBit< BitFEId >, EqBit< BitFEId > >, ordered_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< FiniteElement, boost::string_ref, &FiniteElement::getNameRef > > > > FiniteElement_multiIndex
MultiIndex for entities for FiniteElement.
EntityHandle get_id_for_max_type()
UId getLocalUniqueIdCalculate() const
Generate UId for finite element entity.

◆ getProblemElementsLayout()

MoFEMErrorCode MoFEM::ProblemsManager::getProblemElementsLayout ( const std::string name,
const std::string & fe_name,
PetscLayout * layout ) const

Get layout of elements in the problem.

In layout is stored information how many elements is on each processor, for more information look int petsc documentation http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/IS/PetscLayoutCreate.html#PetscLayoutCreate

Parameters
nameproblem name
fe_namefinite elements
layoutlayout
verbverbosity level
Returns
error code

Definition at line 2593 of file ProblemsManager.cpp.

2595 {
2596 MoFEM::Interface &m_field = cOre;
2597 const Problem *problem_ptr;
2599 CHKERR m_field.get_problem(name, &problem_ptr);
2600 CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2601 fe_name, layout);
2603}
MoFEMErrorCode getNumberOfElementsByNameAndPart(MPI_Comm comm, const std::string name, PetscLayout *layout) const
Get number of finite elements by name on processors.

◆ inheritPartition()

MoFEMErrorCode MoFEM::ProblemsManager::inheritPartition ( const std::string name,
const std::string problem_for_rows,
bool copy_rows,
const std::string problem_for_cols,
bool copy_cols,
int verb = VERBOSE )

build indexing and partition problem inheriting indexing and partitioning from two other problems

Parameters
nameproblem name
problem_for_rowsproblem used to index rows
copy_rowsjust copy rows dofs
problem_for_colsproblem used to index cols
copy_colsjust copy cols dofs

If copy_rows/copy_cols is set to false only partition is copied between problems.

Examples
build_problems.cpp.

Definition at line 1896 of file ProblemsManager.cpp.

1898 {
1899 MoFEM::Interface &m_field = cOre;
1900 auto problems_ptr = m_field.get_problems();
1902
1904 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "pRoblems not build");
1905
1906 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1907
1908 // find p_miit
1909 ProblemByName &problems_by_name =
1910 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1911 ProblemByName::iterator p_miit = problems_by_name.find(name);
1912 if (p_miit == problems_by_name.end()) {
1913 SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
1914 "problem with name < %s > not defined (top tip check spelling)",
1915 name.c_str());
1916 }
1917 if (verb > QUIET)
1918 MOFEM_LOG("WORLD", Sev::inform)
1919 << p_miit->getName() << " from rows of " << problem_for_rows
1920 << " and columns of " << problem_for_cols;
1921
1922 // find p_miit_row
1923 ProblemByName::iterator p_miit_row = problems_by_name.find(problem_for_rows);
1924 if (p_miit_row == problems_by_name.end()) {
1925 SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1926 "problem with name < %s > not defined (top tip check spelling)",
1927 problem_for_rows.c_str());
1928 }
1929 const boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_row =
1930 p_miit_row->numeredRowDofsPtr;
1931
1932 // find p_mit_col
1933 ProblemByName::iterator p_miit_col = problems_by_name.find(problem_for_cols);
1934 if (p_miit_col == problems_by_name.end()) {
1935 SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1936 "problem with name < %s > not defined (top tip check spelling)",
1937 problem_for_cols.c_str());
1938 }
1939 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_col =
1940 p_miit_col->numeredColDofsPtr;
1941
1942 bool copy[] = {copy_rows, copy_cols};
1943 boost::shared_ptr<NumeredDofEntity_multiIndex> composed_dofs[] = {
1944 p_miit->numeredRowDofsPtr, p_miit->numeredColDofsPtr};
1945
1946 int *nb_local_dofs[] = {&p_miit->nbLocDofsRow, &p_miit->nbLocDofsCol};
1947 int *nb_dofs[] = {&p_miit->nbDofsRow, &p_miit->nbDofsCol};
1948 boost::shared_ptr<NumeredDofEntity_multiIndex> copied_dofs[] = {dofs_row,
1949 dofs_col};
1950
1951 for (int ss = 0; ss < 2; ss++) {
1952
1953 // build indices
1954 *nb_local_dofs[ss] = 0;
1955 if (!copy[ss]) {
1956
1957 // only copy indices which are belong to some elements if this problem
1958 std::vector<int> is_local, is_new;
1959
1960 NumeredDofEntityByUId &dofs_by_uid =
1961 copied_dofs[ss]->get<Unique_mi_tag>();
1962 for (NumeredDofEntity_multiIndex::iterator dit =
1963 composed_dofs[ss]->begin();
1964 dit != composed_dofs[ss]->end(); dit++) {
1965 NumeredDofEntityByUId::iterator diit =
1966 dofs_by_uid.find((*dit)->getLocalUniqueId());
1967 if (diit == dofs_by_uid.end()) {
1968 SETERRQ(
1970 "data inconsistency, could not find dof in composite problem");
1971 }
1972 int part_number = (*diit)->getPart(); // get part number
1973 int petsc_global_dof = (*diit)->getPetscGlobalDofIdx();
1974 bool success;
1975 success = composed_dofs[ss]->modify(
1976 dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
1977 petsc_global_dof));
1978 if (!success) {
1979 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1980 "modification unsuccessful");
1981 }
1982 if ((*dit)->getPart() == (unsigned int)m_field.get_comm_rank()) {
1983 success = composed_dofs[ss]->modify(
1984 dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1985 if (!success) {
1986 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1987 "modification unsuccessful");
1988 }
1989 is_local.push_back(petsc_global_dof);
1990 }
1991 }
1992
1993 AO ao;
1994 CHKERR AOCreateMapping(m_field.get_comm(), is_local.size(), &is_local[0],
1995 NULL, &ao);
1996
1997 // apply local to global mapping
1998 is_local.resize(0);
1999 for (NumeredDofEntity_multiIndex::iterator dit =
2000 composed_dofs[ss]->begin();
2001 dit != composed_dofs[ss]->end(); dit++) {
2002 is_local.push_back((*dit)->getPetscGlobalDofIdx());
2003 }
2004 CHKERR AOPetscToApplication(ao, is_local.size(), &is_local[0]);
2005 int idx2 = 0;
2006 for (NumeredDofEntity_multiIndex::iterator dit =
2007 composed_dofs[ss]->begin();
2008 dit != composed_dofs[ss]->end(); dit++) {
2009 int part_number = (*dit)->getPart(); // get part number
2010 int petsc_global_dof = is_local[idx2++];
2011 bool success;
2012 success = composed_dofs[ss]->modify(
2013 dit, NumeredDofEntity_part_and_glob_idx_change(part_number,
2014 petsc_global_dof));
2015 if (!success) {
2016 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2017 "modification unsuccessful");
2018 }
2019 }
2020
2021 CHKERR AODestroy(&ao);
2022
2023 } else {
2024
2025 for (NumeredDofEntity_multiIndex::iterator dit = copied_dofs[ss]->begin();
2026 dit != copied_dofs[ss]->end(); dit++) {
2027 std::pair<NumeredDofEntity_multiIndex::iterator, bool> p;
2028 p = composed_dofs[ss]->insert(boost::shared_ptr<NumeredDofEntity>(
2029 new NumeredDofEntity((*dit)->getDofEntityPtr())));
2030 if (p.second) {
2031 (*nb_dofs[ss])++;
2032 }
2033 int dof_idx = (*dit)->getDofIdx();
2034 int part_number = (*dit)->getPart(); // get part number
2035 int petsc_global_dof = (*dit)->getPetscGlobalDofIdx();
2036 if (part_number == (unsigned int)m_field.get_comm_rank()) {
2037 const bool success = composed_dofs[ss]->modify(
2038 p.first, NumeredDofEntity_part_and_all_indices_change(
2039 part_number, dof_idx, petsc_global_dof,
2040 (*nb_local_dofs[ss])++));
2041 if (!success) {
2042 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2043 "modification unsuccessful");
2044 }
2045 } else {
2046 const bool success = composed_dofs[ss]->modify(
2047 p.first, NumeredDofEntity_part_and_mofem_glob_idx_change(
2048 part_number, dof_idx, petsc_global_dof));
2049 if (!success) {
2050 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2051 "modification unsuccessful");
2052 }
2053 }
2054 }
2055 }
2056 }
2057
2058 CHKERR printPartitionedProblem(&*p_miit, verb);
2059 CHKERR debugPartitionedProblem(&*p_miit, verb);
2060
2062}
@ QUIET
NumeredDofEntity_multiIndex::index< Unique_mi_tag >::type NumeredDofEntityByUId
Numbered DoF multi-index by UId.

◆ partitionFiniteElements()

MoFEMErrorCode MoFEM::ProblemsManager::partitionFiniteElements ( const std::string name,
bool part_from_moab = false,
int low_proc = -1,
int hi_proc = -1,
int verb = VERBOSE )

partition finite elements

Function which partition finite elements based on dofs partitioning.
In addition it sets information about local row and cols dofs at given element on partition.

Parameters
name
part_from_moab
low_proc
hi_proc
verb
Returns
MoFEMErrorCode
Parameters
nameproblem name
Examples
build_problems.cpp, continuity_check_on_skeleton_3d.cpp, forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, hcurl_divergence_operator_2d.cpp, mesh_insert_interface_atom.cpp, nonlinear_dynamics.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, and remove_entities_from_problem.cpp.

Definition at line 2167 of file ProblemsManager.cpp.

2170 {
2171 MoFEM::Interface &m_field = cOre;
2172 auto problems_ptr = m_field.get_problems();
2173 auto fe_ent_ptr = m_field.get_ents_finite_elements();
2175
2177 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "fields not build");
2178
2180 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "FEs not build");
2181
2183 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2184 "adjacencies not build");
2185
2187 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "problem not build");
2188
2190 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2191 "problem not partitioned");
2192
2193 if (low_proc == -1)
2194 low_proc = m_field.get_comm_rank();
2195 if (hi_proc == -1)
2196 hi_proc = m_field.get_comm_rank();
2197
2198 // Find pointer to problem of given name
2199 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
2200 auto &problems =
2201 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
2202 ProblemByName::iterator p_miit = problems.find(name);
2203 if (p_miit == problems.end()) {
2204 SETERRQ1(m_field.get_comm(), MOFEM_NOT_FOUND,
2205 "problem < %s > not found (top tip: check spelling)",
2206 name.c_str());
2207 }
2208
2209 // Get reference on finite elements multi-index on the problem
2210 NumeredEntFiniteElement_multiIndex &problem_finite_elements =
2211 *p_miit->numeredFiniteElementsPtr;
2212
2213 // Clear all elements and data, build it again
2214 problem_finite_elements.clear();
2215
2216 // Check if dofs and columns are the same, i.e. structurally symmetric
2217 // problem
2218 bool do_cols_prob = true;
2219 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
2220 do_cols_prob = false;
2221 }
2222
2223 auto get_good_elems = [&]() {
2224 auto good_elems = std::vector<decltype(fe_ent_ptr->begin())>();
2225 good_elems.reserve(fe_ent_ptr->size());
2226
2227 const auto prb_bit = p_miit->getBitRefLevel();
2228 const auto prb_mask = p_miit->getBitRefLevelMask();
2229
2230 // Loop over all elements in database and if right element is there add it
2231 // to problem finite element multi-index
2232 for (auto efit = fe_ent_ptr->begin(); efit != fe_ent_ptr->end(); ++efit) {
2233
2234 // if element is not part of problem
2235 if (((*efit)->getId() & p_miit->getBitFEId()).any()) {
2236
2237 const auto &fe_bit = (*efit)->getBitRefLevel();
2238
2239 // if entity is not problem refinement level
2240 if ((fe_bit & prb_mask) == fe_bit && (fe_bit & prb_bit).any())
2241 good_elems.emplace_back(efit);
2242 }
2243 }
2244
2245 return good_elems;
2246 };
2247
2248 auto good_elems = get_good_elems();
2249
2250 auto numbered_good_elems_ptr =
2251 boost::make_shared<std::vector<NumeredEntFiniteElement>>();
2252 numbered_good_elems_ptr->reserve(good_elems.size());
2253 for (auto &efit : good_elems)
2254 numbered_good_elems_ptr->emplace_back(NumeredEntFiniteElement(*efit));
2255
2256 if (!do_cols_prob) {
2257 for (auto &fe : *numbered_good_elems_ptr) {
2258 if (fe.sPtr->getRowFieldEntsPtr() == fe.sPtr->getColFieldEntsPtr()) {
2259 fe.getColFieldEntsPtr() = fe.getRowFieldEntsPtr();
2260 }
2261 }
2262 }
2263
2264 if (part_from_moab) {
2265 for (auto &fe : *numbered_good_elems_ptr) {
2266 // if partition is taken from moab partition
2267 int proc = fe.getPartProc();
2268 if (proc == -1 && fe.getEntType() == MBVERTEX)
2269 proc = fe.getOwnerProc();
2270 fe.part = proc;
2271 }
2272 }
2273
2274 for (auto &fe : *numbered_good_elems_ptr) {
2275
2277 CHKERR fe.sPtr->getRowDofView(*(p_miit->numeredRowDofsPtr), rows_view);
2278
2279 if (!part_from_moab) {
2280 std::vector<int> parts(m_field.get_comm_size(), 0);
2281 for (auto &dof_ptr : rows_view)
2282 parts[dof_ptr->pArt]++;
2283 std::vector<int>::iterator pos = max_element(parts.begin(), parts.end());
2284 const auto max_part = std::distance(parts.begin(), pos);
2285 fe.part = max_part;
2286 }
2287 }
2288
2289 for (auto &fe : *numbered_good_elems_ptr) {
2290
2291 auto check_fields_and_dofs = [&]() {
2292 if (!part_from_moab) {
2293 if (fe.getBitFieldIdRow().none() && m_field.get_comm_size() == 0) {
2294 MOFEM_LOG("WORLD", Sev::warning)
2295 << "At least one field has to be added to element row to "
2296 "determine partition of finite element. Check element " +
2297 boost::lexical_cast<std::string>(fe.getName());
2298 }
2299 }
2300
2301 return true;
2302 };
2303
2304 if (check_fields_and_dofs()) {
2305 // Add element to the problem
2306 auto p = problem_finite_elements.insert(
2307 boost::shared_ptr<NumeredEntFiniteElement>(numbered_good_elems_ptr,
2308 &fe));
2309 if (!p.second)
2310 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "element is there");
2311 }
2312 }
2313
2314 if (verb >= VERBOSE) {
2315 auto elements_on_rank =
2316 problem_finite_elements.get<Part_mi_tag>().equal_range(
2317 m_field.get_comm_rank());
2318 MOFEM_LOG("SYNC", Sev::verbose)
2319 << p_miit->getName() << " nb. elems "
2320 << std::distance(elements_on_rank.first, elements_on_rank.second);
2321 auto fe_ptr = m_field.get_finite_elements();
2322 for (auto &fe : *fe_ptr) {
2323 auto e_range =
2324 problem_finite_elements.get<Composite_Name_And_Part_mi_tag>()
2325 .equal_range(
2326 boost::make_tuple(fe->getName(), m_field.get_comm_rank()));
2327 MOFEM_LOG("SYNC", Sev::noisy)
2328 << "Element " << fe->getName() << " nb. elems "
2329 << std::distance(e_range.first, e_range.second);
2330 }
2331
2333 }
2334
2337}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
multi_index_container< boost::shared_ptr< NumeredEntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_EntFiniteElement, UId, &NumeredEntFiniteElement::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > > > > > NumeredEntFiniteElement_multiIndex
MultiIndex for entities for NumeredEntFiniteElement.
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents finite elements object.
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered

◆ partitionGhostDofs()

MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofs ( const std::string name,
int verb = VERBOSE )

determine ghost nodes

Parameters
nameproblem name

DOFs are ghost dofs if are used by elements on given partition, but not owned by that partition.

Examples
build_problems.cpp, continuity_check_on_skeleton_3d.cpp, forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, hcurl_divergence_operator_2d.cpp, mesh_insert_interface_atom.cpp, nonlinear_dynamics.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, and remove_entities_from_problem.cpp.

Definition at line 2339 of file ProblemsManager.cpp.

2340 {
2341 MoFEM::Interface &m_field = cOre;
2342 auto problems_ptr = m_field.get_problems();
2344
2346 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2347 "partition of problem not build");
2349 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2350 "partitions finite elements not build");
2351
2352 // get problem pointer
2353 auto p_miit = problems_ptr->get<Problem_mi_tag>().find(name);
2354 if (p_miit == problems_ptr->get<Problem_mi_tag>().end())
2355 SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Problem %s not fond",
2356 name.c_str());
2357
2358 // get reference to number of ghost dofs
2359 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
2360 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
2361 nb_row_ghost_dofs = 0;
2362 nb_col_ghost_dofs = 0;
2363
2364 // do work if more than one processor
2365 if (m_field.get_comm_size() > 1) {
2366
2368 ghost_idx_col_view;
2369
2370 // get elements on this partition
2371 auto fe_range =
2372 p_miit->numeredFiniteElementsPtr->get<Part_mi_tag>().equal_range(
2373 m_field.get_comm_rank());
2374
2375 // get dofs on elements which are not part of this partition
2376
2377 struct Inserter {
2378 using Vec = std::vector<boost::shared_ptr<NumeredDofEntity>>;
2379 using It = Vec::iterator;
2380 It operator()(Vec &dofs_view, It &hint,
2381 boost::shared_ptr<NumeredDofEntity> &&dof) {
2382 dofs_view.emplace_back(dof);
2383 return dofs_view.end();
2384 }
2385 };
2386
2387 // rows
2388 std::vector<boost::shared_ptr<NumeredDofEntity>> fe_vec_view;
2389 auto hint_r = ghost_idx_row_view.begin();
2390 for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2391
2392 fe_vec_view.clear();
2393 CHKERR EntFiniteElement::getDofView((*fe_ptr)->getRowFieldEnts(),
2394 *(p_miit->getNumeredRowDofsPtr()),
2395 fe_vec_view, Inserter());
2396
2397 for (auto &dof_ptr : fe_vec_view) {
2398 if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2399 hint_r = ghost_idx_row_view.emplace_hint(hint_r, dof_ptr);
2400 }
2401 }
2402 }
2403
2404 // columns
2405 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2406
2407 auto hint_c = ghost_idx_col_view.begin();
2408 for (auto fe_ptr = fe_range.first; fe_ptr != fe_range.second; ++fe_ptr) {
2409
2410 fe_vec_view.clear();
2411 CHKERR EntFiniteElement::getDofView((*fe_ptr)->getColFieldEnts(),
2412 *(p_miit->getNumeredColDofsPtr()),
2413 fe_vec_view, Inserter());
2414
2415 for (auto &dof_ptr : fe_vec_view) {
2416 if (dof_ptr->getPart() != (unsigned int)m_field.get_comm_rank()) {
2417 hint_c = ghost_idx_col_view.emplace_hint(hint_c, dof_ptr);
2418 }
2419 }
2420 }
2421 }
2422
2423 int *nb_ghost_dofs[2] = {&nb_row_ghost_dofs, &nb_col_ghost_dofs};
2424 int nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2425
2427 &ghost_idx_row_view, &ghost_idx_col_view};
2428 NumeredDofEntityByUId *dof_by_uid_no_const[2] = {
2429 &p_miit->numeredRowDofsPtr->get<Unique_mi_tag>(),
2430 &p_miit->numeredColDofsPtr->get<Unique_mi_tag>()};
2431
2432 int loop_size = 2;
2433 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2434 loop_size = 1;
2435 }
2436
2437 // set local ghost dofs indices
2438 for (int ss = 0; ss != loop_size; ++ss) {
2439 for (auto &gid : *ghost_idx_view[ss]) {
2440 NumeredDofEntityByUId::iterator dof =
2441 dof_by_uid_no_const[ss]->find(gid->getLocalUniqueId());
2442 if (PetscUnlikely((*dof)->petscLocalDofIdx != (DofIdx)-1))
2443 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2444 "inconsistent data, ghost dof already set");
2445 bool success = dof_by_uid_no_const[ss]->modify(
2446 dof, NumeredDofEntity_local_idx_change(nb_local_dofs[ss]++));
2447 if (PetscUnlikely(!success))
2448 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
2449 "modification unsuccessful");
2450 (*nb_ghost_dofs[ss])++;
2451 }
2452 }
2453 if (loop_size == 1) {
2454 (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2455 }
2456 }
2457
2458 if (verb > QUIET) {
2459 MOFEM_LOG("SYNC", Sev::inform)
2460 << " FEs ghost dofs on problem " << p_miit->getName()
2461 << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
2462 << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2463 << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
2464
2466 }
2467
2470}
@ MOFEM_INVALID_DATA
Definition definitions.h:36
int DofIdx
Index of DOF.
Definition Types.hpp:18
static MoFEMErrorCode getDofView(const FE_ENTS &fe_ents_view, const MOFEM_DOFS &mofem_dofs, MOFEM_DOFS_VIEW &dofs_view, INSERTER &&inserter)

◆ partitionGhostDofsOnDistributedMesh()

MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofsOnDistributedMesh ( const std::string name,
int verb = VERBOSE )

determine ghost nodes on distributed meshes

Parameters
nameproblem name

It is very similar for partitionGhostDofs, however this explits that mesh is distributed.

DOFs are ghosted if are shered but not owned by partition.

Examples
remove_entities_from_problem.cpp.

Definition at line 2473 of file ProblemsManager.cpp.

2474 {
2475 MoFEM::Interface &m_field = cOre;
2476 auto problems_ptr = m_field.get_problems();
2478
2480 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2481 "partition of problem not build");
2483 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2484 "partitions finite elements not build");
2485
2486 // get problem pointer
2487 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemsByName;
2488 // find p_miit
2489 ProblemsByName &problems_set =
2490 const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
2491 ProblemsByName::iterator p_miit = problems_set.find(name);
2492
2493 // get reference to number of ghost dofs
2494 // get number of local dofs
2495 DofIdx *nb_ghost_dofs[2] = {&(p_miit->nbGhostDofsRow),
2496 &(p_miit->nbGhostDofsCol)};
2497 DofIdx nb_local_dofs[2] = {p_miit->nbLocDofsRow, p_miit->nbLocDofsCol};
2498 for (int ss = 0; ss != 2; ++ss) {
2499 (*nb_ghost_dofs[ss]) = 0;
2500 }
2501
2502 // do work if more than one processor
2503 if (m_field.get_comm_size() > 1) {
2504 // determine if rows on columns are different from dofs on rows
2505 int loop_size = 2;
2506 if (p_miit->numeredColDofsPtr == p_miit->numeredRowDofsPtr) {
2507 loop_size = 1;
2508 }
2509
2510 typedef decltype(p_miit->numeredRowDofsPtr) NumbDofTypeSharedPtr;
2511 NumbDofTypeSharedPtr numered_dofs[] = {p_miit->numeredRowDofsPtr,
2512 p_miit->numeredColDofsPtr};
2513
2514 // iterate over dofs on rows and dofs on columns
2515 for (int ss = 0; ss != loop_size; ++ss) {
2516
2517 // create dofs view by uid
2518 auto r = numered_dofs[ss]->get<PetscLocalIdx_mi_tag>().equal_range(-1);
2519
2520 std::vector<NumeredDofEntity_multiIndex::iterator> ghost_idx_view;
2521 ghost_idx_view.reserve(std::distance(r.first, r.second));
2522 for (; r.first != r.second; ++r.first)
2523 ghost_idx_view.emplace_back(numered_dofs[ss]->project<0>(r.first));
2524
2525 auto cmp = [](auto a, auto b) {
2526 return (*a)->getLocalUniqueId() < (*b)->getLocalUniqueId();
2527 };
2528 sort(ghost_idx_view.begin(), ghost_idx_view.end(), cmp);
2529
2530 // iterate over dofs which have negative local index
2531 for (auto gid_it : ghost_idx_view) {
2532 bool success = numered_dofs[ss]->modify(
2533 gid_it, NumeredDofEntity_local_idx_change((nb_local_dofs[ss])++));
2534 if (!success)
2535 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2536 "modification unsuccessful");
2537 ++(*nb_ghost_dofs[ss]);
2538 }
2539 }
2540 if (loop_size == 1) {
2541 (*nb_ghost_dofs[1]) = (*nb_ghost_dofs[0]);
2542 }
2543 }
2544
2545 if (verb > QUIET) {
2546 MOFEM_LOG("SYNC", Sev::inform)
2547 << " FEs ghost dofs on problem " << p_miit->getName()
2548 << " Nb. ghost dof " << p_miit->getNbGhostDofsRow() << " by "
2549 << p_miit->getNbGhostDofsCol() << " Nb. local dof "
2550 << p_miit->getNbLocalDofsCol() << " by " << p_miit->getNbLocalDofsCol();
2551
2553 }
2554
2557}
constexpr double a

◆ partitionMesh() [1/2]

MoFEMErrorCode MoFEM::CommInterface::partitionMesh ( const Range & ents,
const int dim,
const int adj_dim,
const int n_parts,
Tag * th_vertex_weights = nullptr,
Tag * th_edge_weights = nullptr,
Tag * th_part_weights = nullptr,
int verb = VERBOSE,
const bool debug = false )

Set partition tag to each finite element in the problem.

This will use one of the mesh partitioning programs available from PETSc See http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPartitioningType.html

Parameters
entsEntities to partition
dimDimension of entities to partition
adj_dimAdjacency dimension
n_partsNumber of partitions
verbVerbosity level
Returns
Error code
Examples
partition_mesh.cpp.

Definition at line 875 of file CommInterface.cpp.

878 {
879 MoFEM::Interface &m_field = cOre;
881
882 // get layout
883 int rstart, rend, nb_elems;
884 {
885 PetscLayout layout;
886 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
887 CHKERR PetscLayoutSetBlockSize(layout, 1);
888 CHKERR PetscLayoutSetSize(layout, ents.size());
889 CHKERR PetscLayoutSetUp(layout);
890 CHKERR PetscLayoutGetSize(layout, &nb_elems);
891 CHKERR PetscLayoutGetRange(layout, &rstart, &rend);
892 CHKERR PetscLayoutDestroy(&layout);
893 if (verb >= VERBOSE) {
894 MOFEM_LOG("SYNC", Sev::inform)
895 << "Finite elements in problem: row lower " << rstart << " row upper "
896 << rend << " nb. elems " << nb_elems << " ( " << ents.size() << " )";
898 }
899 }
900
901 std::vector<EntityHandle> weight_ents;
902 weight_ents.reserve(rend - rstart + 1);
903
904 struct AdjBridge {
905 EntityHandle ent;
906 std::vector<int> adj;
907 AdjBridge(const EntityHandle ent, std::vector<int> &adj)
908 : ent(ent), adj(adj) {}
909 };
910
911 typedef multi_index_container<
912 AdjBridge,
913 indexed_by<
914
915 hashed_unique<member<AdjBridge, EntityHandle, &AdjBridge::ent>>
916
917 >>
918 AdjBridgeMap;
919
920 auto get_it = [&](auto i) {
921 auto it = ents.begin();
922 for (; i > 0; --i) {
923 if (it == ents.end())
924 break;
925 ++it;
926 }
927 return it;
928 };
929
930 Range proc_ents;
931 proc_ents.insert(get_it(rstart), get_it(rend));
932 if (proc_ents.size() != rend - rstart)
933 SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
934 "Wrong number of elements in range %d != %d", proc_ents.size(),
935 rend - rstart);
936
937 Range all_dim_ents;
938 CHKERR m_field.get_moab().get_adjacencies(
939 proc_ents, adj_dim, true, all_dim_ents, moab::Interface::UNION);
940
941 AdjBridgeMap adj_bridge_map;
942 auto hint = adj_bridge_map.begin();
943 std::vector<int> adj;
944 for (auto ent : all_dim_ents) {
945 Range adj_ents;
946 CHKERR m_field.get_moab().get_adjacencies(&ent, 1, dim, false, adj_ents);
947 adj_ents = intersect(adj_ents, ents);
948 adj.clear();
949 adj.reserve(adj_ents.size());
950 for (auto a : adj_ents)
951 adj.emplace_back(ents.index(a));
952 hint = adj_bridge_map.emplace_hint(hint, ent, adj);
953 }
954
955 int *_i;
956 int *_j;
957 {
958 const int nb_loc_elements = rend - rstart;
959 std::vector<int> i(nb_loc_elements + 1, 0), j;
960 {
961 std::vector<int> row_adj;
962 Range::iterator fe_it;
963 int ii, jj;
964 size_t max_row_size;
965 for (fe_it = proc_ents.begin(), ii = rstart, jj = 0, max_row_size = 0;
966 fe_it != proc_ents.end(); ++fe_it, ++ii) {
967
968 if (type_from_handle(*fe_it) == MBENTITYSET) {
969 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
970 "not yet implemented, don't know what to do for meshset "
971 "element");
972 } else {
973
974 Range dim_ents;
975 CHKERR m_field.get_moab().get_adjacencies(&*fe_it, 1, adj_dim, false,
976 dim_ents);
977 dim_ents = intersect(dim_ents, all_dim_ents);
978
979 row_adj.clear();
980 for (auto e : dim_ents) {
981 auto adj_it = adj_bridge_map.find(e);
982 if (adj_it != adj_bridge_map.end()) {
983
984 for (const auto idx : adj_it->adj)
985 row_adj.push_back(idx);
986
987 } else
988 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
989 "Entity not found");
990 }
991
992 std::sort(row_adj.begin(), row_adj.end());
993 auto end = std::unique(row_adj.begin(), row_adj.end());
994
995 size_t row_size = std::distance(row_adj.begin(), end);
996 max_row_size = std::max(max_row_size, row_size);
997 if (j.capacity() < (nb_loc_elements - jj) * max_row_size)
998 j.reserve(nb_loc_elements * max_row_size);
999
1000 i[jj] = j.size();
1001 auto diag = ents.index(*fe_it);
1002 for (auto it = row_adj.begin(); it != end; ++it)
1003 if (*it != diag)
1004 j.push_back(*it);
1005 }
1006
1007 ++jj;
1008
1009 if (th_vertex_weights != NULL)
1010 weight_ents.push_back(*fe_it);
1011 }
1012
1013 i[jj] = j.size();
1014 }
1015
1016 CHKERR PetscMalloc(i.size() * sizeof(int), &_i);
1017 CHKERR PetscMalloc(j.size() * sizeof(int), &_j);
1018 copy(i.begin(), i.end(), _i);
1019 copy(j.begin(), j.end(), _j);
1020 }
1021
1022 // get weights
1023 int *vertex_weights = NULL;
1024 if (th_vertex_weights != NULL) {
1025 CHKERR PetscMalloc(weight_ents.size() * sizeof(int), &vertex_weights);
1026 CHKERR m_field.get_moab().tag_get_data(*th_vertex_weights,
1027 &*weight_ents.begin(),
1028 weight_ents.size(), vertex_weights);
1029 }
1030
1031 {
1032 Mat Adj;
1033 // Adjacency matrix used to partition problems, f.e. METIS
1034 CHKERR MatCreateMPIAdj(m_field.get_comm(), rend - rstart, nb_elems, _i, _j,
1035 PETSC_NULL, &Adj);
1036 CHKERR MatSetOption(Adj, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
1037
1038 if (debug) {
1039 Mat A;
1040 CHKERR MatConvert(Adj, MATMPIAIJ, MAT_INITIAL_MATRIX, &A);
1041 CHKERR MatView(A, PETSC_VIEWER_DRAW_WORLD);
1042 std::string wait;
1043 std::cin >> wait;
1044 CHKERR MatDestroy(&A);
1045 }
1046
1047 // run pets to do partitioning
1048 MOFEM_LOG("WORLD", Sev::verbose) << "Start";
1049
1050 MatPartitioning part;
1051 IS is;
1052 CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
1053
1054 CHKERR MatPartitioningSetAdjacency(part, Adj);
1055 CHKERR MatPartitioningSetFromOptions(part);
1056 CHKERR MatPartitioningSetNParts(part, n_parts);
1057 if (th_vertex_weights != NULL) {
1058 CHKERR MatPartitioningSetVertexWeights(part, vertex_weights);
1059 }
1060 PetscBool same;
1061 PetscObjectTypeCompare((PetscObject)part, MATPARTITIONINGPARMETIS, &same);
1062 if (same) {
1063#ifdef PARMETIS
1064 CHKERR MatPartitioningApply_Parmetis_MoFEM(part, &is);
1065#endif
1066 } else {
1067 CHKERR MatPartitioningApply(part, &is);
1068 }
1069
1070 MOFEM_LOG("WORLD", Sev::verbose) << "End";
1071
1072 // gather
1073 IS is_gather, is_num, is_gather_num;
1074 CHKERR ISAllGather(is, &is_gather);
1075 CHKERR ISPartitioningToNumbering(is, &is_num);
1076 CHKERR ISAllGather(is_num, &is_gather_num);
1077
1078 const int *part_number, *gids;
1079 CHKERR ISGetIndices(is_gather, &part_number);
1080 CHKERR ISGetIndices(is_gather_num, &gids);
1081
1082 // set partition tag and gid tag to entities
1083 ParallelComm *pcomm = ParallelComm::get_pcomm(
1084 &m_field.get_moab(), m_field.get_basic_entity_data_ptr()->pcommID);
1085 Tag part_tag = pcomm->part_tag();
1086 CHKERR m_field.get_moab().tag_set_data(part_tag, ents, part_number);
1087 Tag gid_tag = m_field.get_moab().globalId_tag();
1088
1089 std::map<int, Range> parts_ents;
1090 {
1091 // get entities on each part
1092 Range::iterator eit = ents.begin();
1093 for (int ii = 0; eit != ents.end(); eit++, ii++) {
1094 parts_ents[part_number[ii]].insert(*eit);
1095 }
1096 Range tagged_sets;
1097 CHKERR m_field.get_moab().get_entities_by_type_and_tag(
1098 0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
1099 moab::Interface::UNION);
1100 if (!tagged_sets.empty())
1101 CHKERR m_field.get_moab().tag_delete_data(part_tag, tagged_sets);
1102
1103 if (n_parts > (int)tagged_sets.size()) {
1104 // too few partition sets - create missing ones
1105 int num_new = n_parts - tagged_sets.size();
1106 for (int i = 0; i < num_new; i++) {
1107 EntityHandle new_set;
1108 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, new_set);
1109 tagged_sets.insert(new_set);
1110 }
1111 } else if (n_parts < (int)tagged_sets.size()) {
1112 // too many partition sets - delete extras
1113 int num_del = tagged_sets.size() - n_parts;
1114 for (int i = 0; i < num_del; i++) {
1115 EntityHandle old_set = tagged_sets.pop_back();
1116 CHKERR m_field.get_moab().delete_entities(&old_set, 1);
1117 }
1118 }
1119 // write a tag to those sets denoting they're partition sets, with a
1120 // value of the proc number
1121 std::vector<int> dum_ids(n_parts);
1122 for (int i = 0; i < n_parts; i++)
1123 dum_ids[i] = i;
1124 CHKERR m_field.get_moab().tag_set_data(part_tag, tagged_sets,
1125 &*dum_ids.begin());
1126 CHKERR m_field.get_moab().clear_meshset(tagged_sets);
1127
1128 // get lower dimension entities on each part
1129 for (int pp = 0; pp != n_parts; pp++) {
1130 Range dim_ents = parts_ents[pp].subset_by_dimension(dim);
1131 for (int dd = dim - 1; dd >= 0; dd--) {
1132 Range adj_ents;
1133 CHKERR m_field.get_moab().get_adjacencies(
1134 dim_ents, dd, false, adj_ents, moab::Interface::UNION);
1135 parts_ents[pp].merge(adj_ents);
1136 }
1137 }
1138 for (int pp = 1; pp != n_parts; pp++) {
1139 for (int ppp = 0; ppp != pp; ppp++) {
1140 parts_ents[pp] = subtract(parts_ents[pp], parts_ents[ppp]);
1141 }
1142 }
1143
1144 for (int pp = 0; pp != n_parts; pp++) {
1145 CHKERR m_field.get_moab().add_entities(tagged_sets[pp], parts_ents[pp]);
1146 }
1147
1148 auto set_part = [&]() {
1150 for (EntityType t = MBEDGE; t != MBENTITYSET; ++t) {
1151 for (int pp = 0; pp != n_parts; pp++) {
1152 Range type_ents = parts_ents[pp].subset_by_type(t);
1153 CHKERR m_field.get_moab().tag_clear_data(part_tag, type_ents, &pp);
1154 }
1155 }
1157 };
1158
1159 auto set_gid = [&]() {
1161 for (EntityType t = MBVERTEX; t != MBENTITYSET; ++t) {
1162
1163 void *ptr;
1164 int count;
1165
1166 int gid = 1; // moab indexing from 1a
1167 for (int pp = 0; pp != n_parts; pp++) {
1168 Range type_ents = parts_ents[pp].subset_by_type(t);
1169
1170 auto eit = type_ents.begin();
1171 for (; eit != type_ents.end();) {
1172 CHKERR m_field.get_moab().tag_iterate(
1173 gid_tag, eit, type_ents.end(), count, ptr);
1174 auto gid_tag_ptr = static_cast<int *>(ptr);
1175 for (; count > 0; --count) {
1176 *gid_tag_ptr = gid;
1177 ++eit;
1178 ++gid;
1179 ++gid_tag_ptr;
1180 }
1181 }
1182 }
1183 }
1184
1186 };
1187
1188 CHKERR set_part();
1189 CHKERR set_gid();
1190
1191 }
1192
1193 if (debug) {
1194 if (m_field.get_comm_rank() == 0) {
1195 for (int rr = 0; rr != n_parts; rr++) {
1196 ostringstream ss;
1197 ss << "out_part_" << rr << ".vtk";
1198 MOFEM_LOG("SELF", Sev::inform) << "Save debug part mesh " << ss.str();
1199 EntityHandle meshset;
1200 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
1201 CHKERR m_field.get_moab().add_entities(meshset, parts_ents[rr]);
1202 CHKERR m_field.get_moab().write_file(ss.str().c_str(), "VTK", "",
1203 &meshset, 1);
1204 CHKERR m_field.get_moab().delete_entities(&meshset, 1);
1205 }
1206 }
1207 }
1208
1209 CHKERR ISRestoreIndices(is_gather, &part_number);
1210 CHKERR ISRestoreIndices(is_gather_num, &gids);
1211 CHKERR ISDestroy(&is_num);
1212 CHKERR ISDestroy(&is_gather_num);
1213 CHKERR ISDestroy(&is_gather);
1214 CHKERR ISDestroy(&is);
1215 CHKERR MatPartitioningDestroy(&part);
1216 CHKERR MatDestroy(&Adj);
1217 }
1218
1220}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
static const bool debug
constexpr AssemblyType A
constexpr double t
plate stiffness
Definition plate.cpp:58
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.

◆ partitionMesh() [2/2]

MoFEMErrorCode MoFEM::ProblemsManager::partitionMesh ( const Range & ents,
const int dim,
const int adj_dim,
const int n_parts,
Tag * th_vertex_weights = nullptr,
Tag * th_edge_weights = nullptr,
Tag * th_part_weights = nullptr,
int verb = VERBOSE,
const bool debug = false )

Set partition tag to each finite element in the problem.

Deprecated
Use CommInterface

This will use one of the mesh partitioning programs available from PETSc See http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPartitioningType.html

Parameters
entsEntities to partition
dimDimension of entities to partition
adj_dimAdjacency dimension
n_partsNumber of partitions
verbVerbosity level
Returns
Error code

Definition at line 77 of file ProblemsManager.cpp.

80 {
81 return static_cast<MoFEM::Interface &>(cOre)
83 ->partitionMesh(ents, dim, adj_dim, n_parts, th_vertex_weights,
84 th_edge_weights, th_part_weights, verb, debug);
85}
IFACE getInterface() const
Get interface pointer to pointer of interface.

◆ partitionProblem()

MoFEMErrorCode MoFEM::ProblemsManager::partitionProblem ( const std::string name,
int verb = VERBOSE )

partition problem dofs (collective)

Parameters
nameproblem name
Examples
build_problems.cpp, continuity_check_on_skeleton_3d.cpp, mesh_insert_interface_atom.cpp, nonlinear_dynamics.cpp, and remove_entities_from_problem_not_partitioned.cpp.

Definition at line 1683 of file ProblemsManager.cpp.

1684 {
1685 MoFEM::Interface &m_field = cOre;
1686 auto problems_ptr = m_field.get_problems();
1688
1689 MOFEM_LOG("WORLD", Sev::noisy) << "Partition problem " << name;
1690
1691 using NumeredDofEntitysByIdx =
1692 NumeredDofEntity_multiIndex::index<Idx_mi_tag>::type;
1693 using ProblemsByName = Problem_multiIndex::index<Problem_mi_tag>::type;
1694
1695 // Find problem pointer by name
1696 auto &problems_set =
1697 const_cast<ProblemsByName &>(problems_ptr->get<Problem_mi_tag>());
1698 auto p_miit = problems_set.find(name);
1699 if (p_miit == problems_set.end()) {
1700 SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1701 "problem with name %s not defined (top tip check spelling)",
1702 name.c_str());
1703 }
1704 int nb_dofs_row = p_miit->getNbDofsRow();
1705
1706 if (m_field.get_comm_size() != 1) {
1707
1708 if (!(cOre.getBuildMoFEM() & (1 << 0)))
1709 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1710 if (!(cOre.getBuildMoFEM() & (1 << 1)))
1711 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1712 if (!(cOre.getBuildMoFEM() & (1 << 2)))
1713 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1714 "entFEAdjacencies not build");
1715 if (!(cOre.getBuildMoFEM() & (1 << 3)))
1716 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problems not build");
1717
1718 Mat Adj;
1719 CHKERR m_field.getInterface<MatrixManager>()
1720 ->createMPIAdjWithArrays<Idx_mi_tag>(name, &Adj, verb);
1721
1722 int m, n;
1723 CHKERR MatGetSize(Adj, &m, &n);
1724 if (verb > VERY_VERBOSE)
1725 MatView(Adj, PETSC_VIEWER_STDOUT_WORLD);
1726
1727 // partitioning
1728 MatPartitioning part;
1729 IS is;
1730 CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
1731 CHKERR MatPartitioningSetAdjacency(part, Adj);
1732 CHKERR MatPartitioningSetFromOptions(part);
1733 CHKERR MatPartitioningSetNParts(part, m_field.get_comm_size());
1734 CHKERR MatPartitioningApply(part, &is);
1735 if (verb > VERY_VERBOSE)
1736 ISView(is, PETSC_VIEWER_STDOUT_WORLD);
1737
1738 // gather
1739 IS is_gather, is_num, is_gather_num;
1740 CHKERR ISAllGather(is, &is_gather);
1741 CHKERR ISPartitioningToNumbering(is, &is_num);
1742 CHKERR ISAllGather(is_num, &is_gather_num);
1743 const int *part_number, *petsc_idx;
1744 CHKERR ISGetIndices(is_gather, &part_number);
1745 CHKERR ISGetIndices(is_gather_num, &petsc_idx);
1746 int size_is_num, size_is_gather;
1747 CHKERR ISGetSize(is_gather, &size_is_gather);
1748 if (size_is_gather != (int)nb_dofs_row)
1749 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1750 "data inconsistency %d != %d", size_is_gather, nb_dofs_row);
1751
1752 CHKERR ISGetSize(is_num, &size_is_num);
1753 if (size_is_num != (int)nb_dofs_row)
1754 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1755 "data inconsistency %d != %d", size_is_num, nb_dofs_row);
1756
1757 bool square_matrix = false;
1758 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
1759 square_matrix = true;
1760
1761 // if (!square_matrix) {
1762 // // FIXME: This is for back compatibility, if deprecate interface
1763 // function
1764 // // build interfaces is removed, this part of the code will be obsolete
1765 // auto mit_row = p_miit->numeredRowDofsPtr->get<Idx_mi_tag>().begin();
1766 // auto hi_mit_row = p_miit->numeredRowDofsPtr->get<Idx_mi_tag>().end();
1767 // auto mit_col = p_miit->numeredColDofsPtr->get<Idx_mi_tag>().begin();
1768 // auto hi_mit_col = p_miit->numeredColDofsPtr->get<Idx_mi_tag>().end();
1769 // for (; mit_row != hi_mit_row; mit_row++, mit_col++) {
1770 // if (mit_col == hi_mit_col) {
1771 // SETERRQ(
1772 // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1773 // "check finite element definition, nb. of rows is not equal to "
1774 // "number for columns");
1775 // }
1776 // if (mit_row->get()->getLocalUniqueId() !=
1777 // mit_col->get()->getLocalUniqueId()) {
1778 // SETERRQ(
1779 // PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1780 // "check finite element definition, nb. of rows is not equal to "
1781 // "number for columns");
1782 // }
1783 // }
1784 // }
1785
1786 auto number_dofs = [&](auto &dofs_idx, auto &counter) {
1788 for (auto miit_dofs_row = dofs_idx.begin();
1789 miit_dofs_row != dofs_idx.end(); miit_dofs_row++) {
1790 const int part = part_number[(*miit_dofs_row)->dofIdx];
1791 if (part == (unsigned int)m_field.get_comm_rank()) {
1792 const bool success = dofs_idx.modify(
1793 miit_dofs_row,
1794 NumeredDofEntity_part_and_indices_change(
1795 part, petsc_idx[(*miit_dofs_row)->dofIdx], counter++));
1796 if (!success) {
1797 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1798 "modification unsuccessful");
1799 }
1800 } else {
1801 const bool success = dofs_idx.modify(
1802 miit_dofs_row, NumeredDofEntity_part_and_glob_idx_change(
1803 part, petsc_idx[(*miit_dofs_row)->dofIdx]));
1804 if (!success) {
1805 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1806 "modification unsuccessful");
1807 }
1808 }
1809 }
1811 };
1812
1813 // Set petsc global indices
1814 auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
1815 p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
1816 int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1817 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1818 nb_row_local_dofs = 0;
1819 nb_row_ghost_dofs = 0;
1820
1821 CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
1822
1823 int &nb_col_local_dofs = p_miit->nbLocDofsCol;
1824 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1825 if (square_matrix) {
1826 nb_col_local_dofs = nb_row_local_dofs;
1827 nb_col_ghost_dofs = nb_row_ghost_dofs;
1828 } else {
1829 NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
1830 const_cast<NumeredDofEntitysByIdx &>(
1831 p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
1832 nb_col_local_dofs = 0;
1833 nb_col_ghost_dofs = 0;
1834 CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
1835 }
1836
1837 CHKERR ISRestoreIndices(is_gather, &part_number);
1838 CHKERR ISRestoreIndices(is_gather_num, &petsc_idx);
1839 CHKERR ISDestroy(&is_num);
1840 CHKERR ISDestroy(&is_gather_num);
1841 CHKERR ISDestroy(&is_gather);
1842 CHKERR ISDestroy(&is);
1843 CHKERR MatPartitioningDestroy(&part);
1844 CHKERR MatDestroy(&Adj);
1845 CHKERR printPartitionedProblem(&*p_miit, verb);
1846 } else {
1847
1848 auto number_dofs = [&](auto &dof_idx, auto &counter) {
1850 for (auto miit_dofs_row = dof_idx.begin(); miit_dofs_row != dof_idx.end();
1851 miit_dofs_row++) {
1852 const bool success = dof_idx.modify(
1853 miit_dofs_row,
1854 NumeredDofEntity_part_and_indices_change(0, counter, counter));
1855 ++counter;
1856 if (!success) {
1857 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1858 "modification unsuccessful");
1859 }
1860 }
1862 };
1863
1864 auto &dofs_row_by_idx_no_const = const_cast<NumeredDofEntitysByIdx &>(
1865 p_miit->numeredRowDofsPtr->get<Idx_mi_tag>());
1866 int &nb_row_local_dofs = p_miit->nbLocDofsRow;
1867 int &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1868 nb_row_local_dofs = 0;
1869 nb_row_ghost_dofs = 0;
1870
1871 CHKERR number_dofs(dofs_row_by_idx_no_const, nb_row_local_dofs);
1872
1873 bool square_matrix = false;
1874 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr)
1875 square_matrix = true;
1876
1877 int &nb_col_local_dofs = p_miit->nbLocDofsCol;
1878 int &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1879 if (square_matrix) {
1880 nb_col_local_dofs = nb_row_local_dofs;
1881 nb_col_ghost_dofs = nb_row_ghost_dofs;
1882 } else {
1883 NumeredDofEntitysByIdx &dofs_col_by_idx_no_const =
1884 const_cast<NumeredDofEntitysByIdx &>(
1885 p_miit->numeredColDofsPtr->get<Idx_mi_tag>());
1886 nb_col_local_dofs = 0;
1887 nb_col_ghost_dofs = 0;
1888 CHKERR number_dofs(dofs_col_by_idx_no_const, nb_col_local_dofs);
1889 }
1890 }
1891
1892 cOre.getBuildMoFEM() |= 1 << 4;
1894}
@ VERY_VERBOSE
const double n
refractive index of diffusive medium
FTensor::Index< 'm', 3 > m

◆ partitionSimpleProblem()

MoFEMErrorCode MoFEM::ProblemsManager::partitionSimpleProblem ( const std::string name,
int verb = VERBOSE )

partition problem dofs

Parameters
nameproblem name
Examples
forces_and_sources_testing_edge_element.cpp, forces_and_sources_testing_flat_prism_element.cpp, hcurl_divergence_operator_2d.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, and quad_polynomial_approximation.cpp.

Definition at line 1537 of file ProblemsManager.cpp.

1538 {
1539
1540 MoFEM::Interface &m_field = cOre;
1541 auto problems_ptr = m_field.get_problems();
1543
1545 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
1547 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
1549 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1551 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
1552 MOFEM_LOG("WORLD", Sev::verbose) << "Simple partition problem " << name;
1553
1554 // find p_miit
1555 typedef Problem_multiIndex::index<Problem_mi_tag>::type ProblemByName;
1556 ProblemByName &problems_set =
1557 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1558 ProblemByName::iterator p_miit = problems_set.find(name);
1559 if (p_miit == problems_set.end()) {
1560 SETERRQ1(PETSC_COMM_SELF, 1,
1561 "problem < %s > is not found (top tip: check spelling)",
1562 name.c_str());
1563 }
1564 typedef boost::multi_index::index<NumeredDofEntity_multiIndex,
1565 Idx_mi_tag>::type NumeredDofEntitysByIdx;
1566 NumeredDofEntitysByIdx &dofs_row_by_idx =
1567 p_miit->numeredRowDofsPtr->get<Idx_mi_tag>();
1568 NumeredDofEntitysByIdx &dofs_col_by_idx =
1569 p_miit->numeredColDofsPtr->get<Idx_mi_tag>();
1570 boost::multi_index::index<NumeredDofEntity_multiIndex,
1571 Idx_mi_tag>::type::iterator miit_row,
1572 hi_miit_row;
1573 boost::multi_index::index<NumeredDofEntity_multiIndex,
1574 Idx_mi_tag>::type::iterator miit_col,
1575 hi_miit_col;
1576 DofIdx &nb_row_local_dofs = p_miit->nbLocDofsRow;
1577 DofIdx &nb_row_ghost_dofs = p_miit->nbGhostDofsRow;
1578 nb_row_local_dofs = 0;
1579 nb_row_ghost_dofs = 0;
1580 DofIdx &nb_col_local_dofs = p_miit->nbLocDofsCol;
1581 DofIdx &nb_col_ghost_dofs = p_miit->nbGhostDofsCol;
1582 nb_col_local_dofs = 0;
1583 nb_col_ghost_dofs = 0;
1584
1585 bool square_matrix = false;
1586 if (p_miit->numeredRowDofsPtr == p_miit->numeredColDofsPtr) {
1587 square_matrix = true;
1588 }
1589
1590 // get row range of local indices
1591 PetscLayout layout_row;
1592 const int *ranges_row;
1593
1594 DofIdx nb_dofs_row = dofs_row_by_idx.size();
1595 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_row);
1596 CHKERR PetscLayoutSetBlockSize(layout_row, 1);
1597 CHKERR PetscLayoutSetSize(layout_row, nb_dofs_row);
1598 CHKERR PetscLayoutSetUp(layout_row);
1599 CHKERR PetscLayoutGetRanges(layout_row, &ranges_row);
1600 // get col range of local indices
1601 PetscLayout layout_col;
1602 const int *ranges_col;
1603 if (!square_matrix) {
1604 DofIdx nb_dofs_col = dofs_col_by_idx.size();
1605 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout_col);
1606 CHKERR PetscLayoutSetBlockSize(layout_col, 1);
1607 CHKERR PetscLayoutSetSize(layout_col, nb_dofs_col);
1608 CHKERR PetscLayoutSetUp(layout_col);
1609 CHKERR PetscLayoutGetRanges(layout_col, &ranges_col);
1610 }
1611 for (unsigned int part = 0; part < (unsigned int)m_field.get_comm_size();
1612 part++) {
1613 miit_row = dofs_row_by_idx.lower_bound(ranges_row[part]);
1614 hi_miit_row = dofs_row_by_idx.lower_bound(ranges_row[part + 1]);
1615 if (std::distance(miit_row, hi_miit_row) !=
1616 ranges_row[part + 1] - ranges_row[part]) {
1617 SETERRQ4(
1618 PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1619 "data inconsistency, std::distance(miit_row,hi_miit_row) != rend - "
1620 "rstart (%d != %d - %d = %d) ",
1621 std::distance(miit_row, hi_miit_row), ranges_row[part + 1],
1622 ranges_row[part], ranges_row[part + 1] - ranges_row[part]);
1623 }
1624 // loop rows
1625 for (; miit_row != hi_miit_row; miit_row++) {
1626 bool success = dofs_row_by_idx.modify(
1627 miit_row,
1628 NumeredDofEntity_part_and_glob_idx_change(part, (*miit_row)->dofIdx));
1629 if (!success)
1630 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1631 "modification unsuccessful");
1632 if (part == (unsigned int)m_field.get_comm_rank()) {
1633 success = dofs_row_by_idx.modify(
1634 miit_row, NumeredDofEntity_local_idx_change(nb_row_local_dofs++));
1635 if (!success)
1636 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1637 "modification unsuccessful");
1638 }
1639 }
1640 if (!square_matrix) {
1641 miit_col = dofs_col_by_idx.lower_bound(ranges_col[part]);
1642 hi_miit_col = dofs_col_by_idx.lower_bound(ranges_col[part + 1]);
1643 if (std::distance(miit_col, hi_miit_col) !=
1644 ranges_col[part + 1] - ranges_col[part]) {
1645 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1646 "data inconsistency, std::distance(miit_col,hi_miit_col) != "
1647 "rend - "
1648 "rstart (%d != %d - %d = %d) ",
1649 std::distance(miit_col, hi_miit_col), ranges_col[part + 1],
1650 ranges_col[part], ranges_col[part + 1] - ranges_col[part]);
1651 }
1652 // loop cols
1653 for (; miit_col != hi_miit_col; miit_col++) {
1654 bool success = dofs_col_by_idx.modify(
1655 miit_col, NumeredDofEntity_part_and_glob_idx_change(
1656 part, (*miit_col)->dofIdx));
1657 if (!success)
1658 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1659 "modification unsuccessful");
1660 if (part == (unsigned int)m_field.get_comm_rank()) {
1661 success = dofs_col_by_idx.modify(
1662 miit_col, NumeredDofEntity_local_idx_change(nb_col_local_dofs++));
1663 if (!success)
1664 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1665 "modification unsuccessful");
1666 }
1667 }
1668 }
1669 }
1670 CHKERR PetscLayoutDestroy(&layout_row);
1671 if (!square_matrix) {
1672 CHKERR PetscLayoutDestroy(&layout_col);
1673 }
1674 if (square_matrix) {
1675 nb_col_local_dofs = nb_row_local_dofs;
1676 nb_col_ghost_dofs = nb_row_ghost_dofs;
1677 }
1678 CHKERR printPartitionedProblem(&*p_miit, verb);
1681}

◆ removeDofsOnEntities() [1/2]

MoFEMErrorCode MoFEM::ProblemsManager::removeDofsOnEntities ( const std::string problem_name,
const std::string field_name,
const BitRefLevel bit_ref_level,
const BitRefLevel bit_ref_mask,
Range * ents_ptr = nullptr,
const int lo_coeff = 0,
const int hi_coeff = MAX_DOFS_ON_ENTITY,
const int lo_order = 0,
const int hi_order = 100,
int verb = VERBOSE,
const bool debug = false )

Remove DOFs from problem by bit ref level.

See for more detail other implementation for removeDofsOnEntities.

Parameters
problem_namename of the problem
field_namename of the field
bit_ref_levelbit ref level on which DOFs are removed
bit_ref_maskbit ref mask on which DOFs are removed
ents_ptrfilter entities with given bit and mask
lo_coefflow dof coefficient (rank)
hi_coeffhigh dof coefficient (rank)
verbverbosity level
debugto debug and seek for inconsistencies set to true
Returns
MoFEMErrorCode

Definition at line 3494 of file ProblemsManager.cpp.

3498 {
3499 MoFEM::Interface &m_field = cOre;
3501
3502 auto bit_manager = m_field.getInterface<BitRefManager>();
3503
3504 Range ents;
3505 if (ents_ptr) {
3506 ents = *ents_ptr;
3507 CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3508 ents, verb);
3509 } else {
3510 CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3511 verb);
3512 }
3513
3514 CHKERR removeDofsOnEntities(problem_name, field_name, ents, lo_coeff,
3515 hi_coeff, lo_order, hi_order, verb, debug);
3516
3518}
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
constexpr auto field_name

◆ removeDofsOnEntities() [2/2]

MoFEMErrorCode MoFEM::ProblemsManager::removeDofsOnEntities ( const std::string problem_name,
const std::string field_name,
const Range ents,
const int lo_coeff = 0,
const int hi_coeff = MAX_DOFS_ON_ENTITY,
const int lo_order = 0,
const int hi_order = 100,
int verb = VERBOSE,
const bool debug = false )

Remove DOFs from problem.

Remove DOFs from problem which are on entities on the given range and given field name. On the finite element level, DOFs can be still accessed however local PETSc indices and global PETSc indices are marked with the index -1.

Note
If the index is marked -1 it is not assembled and dropped by VecSetValues and MatSetValues.
Todo
Not yet implemented update for AO maps and IS ranges if removed entities in composite problem or sub-problem
Parameters
problem_namename of the problem
field_namename of the field
entsentities on which DOFs are removed
lo_coefflow dof coefficient (rank)
hi_coeffhigh dof coefficient (rank)
verbverbosity level
debugto debug and seek for inconsistencies set to true
Returns
MoFEMErrorCode
Examples
hanging_node_approx.cpp, helmholtz.cpp, and remove_entities_from_problem.cpp.

Definition at line 2954 of file ProblemsManager.cpp.

2957 {
2958
2959 MoFEM::Interface &m_field = cOre;
2961
2962 const Problem *prb_ptr;
2963 CHKERR m_field.get_problem(problem_name, &prb_ptr);
2964
2965 decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
2966 prb_ptr->numeredRowDofsPtr, nullptr};
2967 if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2968 numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2969
2970 int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2971 int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2972 int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2973
2974 const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
2975 const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
2976 const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
2977 // const int nb_init_loc_col_dofs = prb_ptr->getNbLocalDofsCol();
2978 const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
2979 const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
2980
2981 for (int s = 0; s != 2; ++s)
2982 if (numered_dofs[s]) {
2983
2984 using NumeredDofEntity_it_view_multiIndex = multi_index_container<
2985
2986 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2987
2988 >;
2989
2990 const auto bit_number = m_field.get_field_bit_number(field_name);
2991 NumeredDofEntity_it_view_multiIndex dofs_it_view;
2992
2993 // Set -1 to global and local dofs indices
2994 for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2995 ++pit) {
2996 auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
2997 DofEntity::getLoFieldEntityUId(bit_number, pit->first));
2998 auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
2999 DofEntity::getHiFieldEntityUId(bit_number, pit->second));
3000
3001 for (; lo != hi; ++lo)
3002 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
3003 (*lo)->getDofCoeffIdx() <= hi_coeff &&
3004 (*lo)->getDofOrder() >= lo_order &&
3005 (*lo)->getDofOrder() <= hi_order)
3006 dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
3007 }
3008
3009 if (verb > QUIET) {
3010 for (auto &dof : dofs_it_view)
3011 MOFEM_LOG("SYNC", Sev::noisy) << **dof;
3012 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
3013 }
3014
3015 // create weak view
3016 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_weak_view;
3017 dofs_weak_view.reserve(dofs_it_view.size());
3018 for (auto dit : dofs_it_view)
3019 dofs_weak_view.push_back(*dit);
3020
3021 if (verb >= NOISY)
3022 MOFEM_LOG_C("SYNC", Sev::noisy,
3023 "Number of DOFs in multi-index %d and to delete %d\n",
3024 numered_dofs[s]->size(), dofs_it_view.size());
3025
3026 // erase dofs from problem
3027 for (auto weak_dit : dofs_weak_view)
3028 if (auto dit = weak_dit.lock()) {
3029 numered_dofs[s]->erase(dit->getLocalUniqueId());
3030 }
3031
3032 if (verb >= NOISY)
3033 MOFEM_LOG_C("SYNC", Sev::noisy,
3034 "Number of DOFs in multi-index after delete %d\n",
3035 numered_dofs[s]->size());
3036
3037 // get current number of ghost dofs
3038 int nb_local_dofs = 0;
3039 int nb_ghost_dofs = 0;
3040 for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
3041 dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
3042 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
3043 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
3044 ++nb_local_dofs;
3045 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
3046 ++nb_ghost_dofs;
3047 else
3048 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3049 "Impossible case. You could run problem on no distributed "
3050 "mesh. That is not implemented. Dof local index is %d",
3051 (*dit)->getPetscLocalDofIdx());
3052 }
3053
3054 // get indices
3055 auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
3056 const int nb_dofs = numered_dofs[s]->size();
3057 indices.clear();
3058 indices.reserve(nb_dofs);
3059 for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
3060 dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
3061 bool add = true;
3062 if (only_local) {
3063 if ((*dit)->getPetscLocalDofIdx() < 0 ||
3064 (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
3065 add = false;
3066 }
3067 }
3068 if (add)
3069 indices.push_back(decltype(tag)::get_index(dit));
3070 }
3071 };
3072
3073 auto get_indices_by_uid = [&](auto tag, auto &indices) {
3074 const int nb_dofs = numered_dofs[s]->size();
3075 indices.clear();
3076 indices.reserve(nb_dofs);
3077 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3078 ++dit)
3079 indices.push_back(decltype(tag)::get_index(dit));
3080 };
3081
3082 auto get_sub_ao = [&](auto sub_data) {
3083 if (s == 0) {
3084 return sub_data->getSmartRowMap();
3085 } else {
3086 return sub_data->getSmartColMap();
3087 }
3088 };
3089
3090 auto set_sub_is_and_ao = [&s, &prb_ptr](auto sub_data, auto is, auto ao) {
3091 if (s == 0) {
3092 sub_data->rowIs = is;
3093 sub_data->rowMap = ao;
3094 sub_data->colIs = is; // if square problem col is equal to row
3095 sub_data->colMap = ao;
3096 } else {
3097 sub_data->colIs = is;
3098 sub_data->colMap = ao;
3099 }
3100 };
3101
3102 auto apply_symmetry = [&s, &prb_ptr](auto sub_data) {
3103 if (s == 0) {
3104 if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
3105 sub_data->colIs = sub_data->getSmartRowIs();
3106 sub_data->colMap = sub_data->getSmartRowMap();
3107 }
3108 }
3109 };
3110
3111 auto concatenate_dofs = [&](auto tag, auto &indices,
3112 const auto local_only) {
3114 get_indices_by_tag(tag, indices, local_only);
3115
3116 SmartPetscObj<AO> ao;
3117 // Create AO from app indices (i.e. old), to pestc indices (new after
3118 // remove)
3119 if (local_only)
3120 ao = createAOMapping(m_field.get_comm(), indices.size(),
3121 &*indices.begin(), PETSC_NULL);
3122 else
3123 ao = createAOMapping(PETSC_COMM_SELF, indices.size(),
3124 &*indices.begin(), PETSC_NULL);
3125
3126 // Set mapping to sub dm data
3127 if (local_only) {
3128 if (auto sub_data = prb_ptr->getSubData()) {
3129 // create is and then map it to main problem of sub-problem
3130 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
3131 &*indices.begin(), PETSC_COPY_VALUES);
3132 // get old app, i.e. oroginal befor sub indices, and ao, from app,
3133 // to petsc sub indices.
3134 auto sub_ao = get_sub_ao(sub_data);
3135 CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
3136 sub_ao = createAOMappingIS(sub_is, PETSC_NULL);
3137 // set new sub ao
3138 set_sub_is_and_ao(sub_data, sub_is, sub_ao);
3139 apply_symmetry(sub_data);
3140 } else {
3141 // create sub data
3142 prb_ptr->getSubData() =
3143 boost::make_shared<Problem::SubProblemData>();
3144 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
3145 &*indices.begin(), PETSC_COPY_VALUES);
3146 // set sub is ao
3147 set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, ao);
3148 apply_symmetry(prb_ptr->getSubData());
3149 }
3150 }
3151
3152 get_indices_by_uid(tag, indices);
3153 CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
3154
3156 };
3157
3158 // set indices index
3159 auto set_concatenated_indices = [&]() {
3160 std::vector<int> global_indices;
3161 std::vector<int> local_indices;
3163 CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
3164 CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
3165 auto gi = global_indices.begin();
3166 auto li = local_indices.begin();
3167 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
3168 ++dit) {
3169 auto mod = NumeredDofEntity_part_and_all_indices_change(
3170 (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
3171 bool success = numered_dofs[s]->modify(dit, mod);
3172 if (!success)
3173 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
3174 "can not set negative indices");
3175 ++gi;
3176 ++li;
3177 }
3179 };
3180 CHKERR set_concatenated_indices();
3181
3182 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
3183 m_field.get_comm());
3184 *(local_nbdof_ptr[s]) = nb_local_dofs;
3185 *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
3186
3187 if (debug)
3188 for (auto dof : (*numered_dofs[s])) {
3189 if (dof->getPetscGlobalDofIdx() < 0) {
3190 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3191 "Negative global idx");
3192 }
3193 if (dof->getPetscLocalDofIdx() < 0) {
3194 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
3195 "Negative local idx");
3196 }
3197 }
3198
3199 } else {
3200
3201 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
3202 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
3203 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
3204 }
3205
3206 if (verb > QUIET) {
3208 "WORLD", Sev::inform,
3209 "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
3210 prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
3211 prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
3212 MOFEM_LOG_C("SYNC", Sev::verbose,
3213 "Removed DOFs from problem %s dofs [ %d / %d "
3214 "(before %d / %d) local, %d / %d (before %d / %d)]",
3215 prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
3216 prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
3217 nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
3218 prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
3219 nb_init_ghost_col_dofs);
3220 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
3221 }
3222
3224}
#define MOFEM_LOG_C(channel, severity, format,...)
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.
auto createAOMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[])
Creates an application mapping using two integer arrays.
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredColDofsPtr
store DOFs on columns for this problem