v0.14.0
Loading...
Searching...
No Matches
Files | Classes | Functions
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
MoFEM::LogManager::SeverityLevel Sev
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 refernce 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}
static Index< 'p', 3 > p
@ 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.
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 868 of file CommInterface.cpp.

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

◆ 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
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > 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 3142 of file ProblemsManager.cpp.

3146 {
3147 MoFEM::Interface &m_field = cOre;
3149
3150 auto bit_manager = m_field.getInterface<BitRefManager>();
3151
3152 Range ents;
3153 if (ents_ptr) {
3154 ents = *ents_ptr;
3155 CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3156 ents, verb);
3157 } else {
3158 CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3159 verb);
3160 }
3161
3162 CHKERR removeDofsOnEntities(problem_name, field_name, ents, lo_coeff,
3163 hi_coeff, lo_order, hi_order, verb, debug);
3164
3166}
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
EshelbianPlasticity.cpp, hanging_node_approx.cpp, helmholtz.cpp, and remove_entities_from_problem.cpp.

Definition at line 2605 of file ProblemsManager.cpp.

2608 {
2609
2610 MoFEM::Interface &m_field = cOre;
2612
2613 const Problem *prb_ptr;
2614 CHKERR m_field.get_problem(problem_name, &prb_ptr);
2615
2616 decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
2617 prb_ptr->numeredRowDofsPtr, nullptr};
2618 if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2619 numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2620
2621 int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2622 int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2623 int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2624
2625 const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
2626 const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
2627 const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
2628 // const int nb_init_loc_col_dofs = prb_ptr->getNbLocalDofsCol();
2629 const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
2630 const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
2631
2632 for (int s = 0; s != 2; ++s)
2633 if (numered_dofs[s]) {
2634
2635 typedef multi_index_container<
2636
2637 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2638
2639 >
2640 NumeredDofEntity_it_view_multiIndex;
2641
2642 const auto bit_number = m_field.get_field_bit_number(field_name);
2643 NumeredDofEntity_it_view_multiIndex dofs_it_view;
2644
2645 // Set -1 to global and local dofs indices
2646 for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2647 ++pit) {
2648 auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
2649 DofEntity::getLoFieldEntityUId(bit_number, pit->first));
2650 auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
2651 DofEntity::getHiFieldEntityUId(bit_number, pit->second));
2652
2653 for (; lo != hi; ++lo)
2654 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2655 (*lo)->getDofCoeffIdx() <= hi_coeff &&
2656 (*lo)->getDofOrder() >= lo_order &&
2657 (*lo)->getDofOrder() <= hi_order)
2658 dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2659 }
2660
2661 if (verb > QUIET) {
2662 for (auto &dof : dofs_it_view)
2663 MOFEM_LOG("SYNC", Sev::noisy) << **dof;
2664 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
2665 }
2666
2667 // create weak view
2668 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_weak_view;
2669 dofs_weak_view.reserve(dofs_it_view.size());
2670 for (auto dit : dofs_it_view)
2671 dofs_weak_view.push_back(*dit);
2672
2673 if (verb >= NOISY)
2674 MOFEM_LOG_C("SYNC", Sev::noisy,
2675 "Number of DOFs in multi-index %d and to delete %d\n",
2676 numered_dofs[s]->size(), dofs_it_view.size());
2677
2678 // erase dofs from problem
2679 for (auto weak_dit : dofs_weak_view)
2680 if (auto dit = weak_dit.lock()) {
2681 numered_dofs[s]->erase(dit->getLocalUniqueId());
2682 }
2683
2684 if (verb >= NOISY)
2685 MOFEM_LOG_C("SYNC", Sev::noisy,
2686 "Number of DOFs in multi-index after delete %d\n",
2687 numered_dofs[s]->size());
2688
2689 // get current number of ghost dofs
2690 int nb_local_dofs = 0;
2691 int nb_ghost_dofs = 0;
2692 for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
2693 dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
2694 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2695 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2696 ++nb_local_dofs;
2697 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2698 ++nb_ghost_dofs;
2699 else
2700 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2701 "Impossible case. You could run problem on no distributed "
2702 "mesh. That is not implemented. Dof local index is %d",
2703 (*dit)->getPetscLocalDofIdx());
2704 }
2705
2706 // get indices
2707 auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
2708 const int nb_dofs = numered_dofs[s]->size();
2709 indices.clear();
2710 indices.reserve(nb_dofs);
2711 for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
2712 dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
2713 bool add = true;
2714 if (only_local) {
2715 if ((*dit)->getPetscLocalDofIdx() < 0 ||
2716 (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
2717 add = false;
2718 }
2719 }
2720 if (add)
2721 indices.push_back(decltype(tag)::get_index(dit));
2722 }
2723 };
2724
2725 auto get_indices_by_uid = [&](auto tag, auto &indices) {
2726 const int nb_dofs = numered_dofs[s]->size();
2727 indices.clear();
2728 indices.reserve(nb_dofs);
2729 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2730 ++dit)
2731 indices.push_back(decltype(tag)::get_index(dit));
2732 };
2733
2734 auto get_sub_ao = [&](auto sub_data) {
2735 if (s == 0) {
2736 return sub_data->getSmartRowMap();
2737 } else {
2738 return sub_data->getSmartColMap();
2739 }
2740 };
2741
2742 auto set_sub_is_and_ao = [&s, &prb_ptr](auto sub_data, auto is, auto ao) {
2743 if (s == 0) {
2744 sub_data->rowIs = is;
2745 sub_data->rowMap = ao;
2746 } else {
2747 sub_data->colIs = is;
2748 sub_data->colMap = ao;
2749 }
2750 };
2751
2752 auto apply_symmetry = [&s, &prb_ptr](auto sub_data) {
2753 if (s == 0) {
2754 if (prb_ptr->numeredRowDofsPtr == prb_ptr->numeredColDofsPtr) {
2755 sub_data->colIs = sub_data->getSmartRowIs();
2756 sub_data->colMap = sub_data->getSmartRowMap();
2757 }
2758 }
2759 };
2760
2761 auto concatenate_dofs = [&](auto tag, auto &indices,
2762 const auto local_only) {
2764 get_indices_by_tag(tag, indices, local_only);
2765
2766 SmartPetscObj<AO> ao;
2767 // Create AO from app indices (i.e. old), to pestc indices (new after
2768 // remove)
2769 if (local_only)
2770 ao = createAOMapping(m_field.get_comm(), indices.size(),
2771 &*indices.begin(), PETSC_NULL);
2772 else
2773 ao = createAOMapping(PETSC_COMM_SELF, indices.size(),
2774 &*indices.begin(), PETSC_NULL);
2775
2776 // Set mapping to sub dm data
2777 if (local_only) {
2778 if (auto sub_data = prb_ptr->getSubData()) {
2779 // create is and then map it to main problem of sub-problem
2780 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
2781 &*indices.begin(), PETSC_COPY_VALUES);
2782 // get old app, i.e. oroginal befor sub indices, and ao, from app,
2783 // to petsc sub indices.
2784 auto sub_ao = get_sub_ao(sub_data);
2785 CHKERR AOPetscToApplicationIS(sub_ao, sub_is);
2786 sub_ao = createAOMappingIS(sub_is, PETSC_NULL);
2787 // set new sub ao
2788 set_sub_is_and_ao(sub_data, sub_is, sub_ao);
2789 apply_symmetry(sub_data);
2790 } else {
2791 // create sub data
2792 prb_ptr->getSubData() =
2793 boost::make_shared<Problem::SubProblemData>();
2794 auto sub_is = createISGeneral(m_field.get_comm(), indices.size(),
2795 &*indices.begin(), PETSC_COPY_VALUES);
2796 // set sub is ao
2797 set_sub_is_and_ao(prb_ptr->getSubData(), sub_is, ao);
2798 apply_symmetry(prb_ptr->getSubData());
2799 }
2800 }
2801
2802 get_indices_by_uid(tag, indices);
2803 CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
2804
2806 };
2807
2808 // set indices index
2809 auto set_concatenated_indices = [&]() {
2810 std::vector<int> global_indices;
2811 std::vector<int> local_indices;
2813 CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
2814 CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
2815 auto gi = global_indices.begin();
2816 auto li = local_indices.begin();
2817 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2818 ++dit) {
2819 auto mod = NumeredDofEntity_part_and_all_indices_change(
2820 (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
2821 bool success = numered_dofs[s]->modify(dit, mod);
2822 if (!success)
2823 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2824 "can not set negative indices");
2825 ++gi;
2826 ++li;
2827 }
2829 };
2830 CHKERR set_concatenated_indices();
2831
2832 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
2833 m_field.get_comm());
2834 *(local_nbdof_ptr[s]) = nb_local_dofs;
2835 *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
2836
2837 if (debug)
2838 for (auto dof : (*numered_dofs[s])) {
2839 if (dof->getPetscGlobalDofIdx() < 0) {
2840 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2841 "Negative global idx");
2842 }
2843 if (dof->getPetscLocalDofIdx() < 0) {
2844 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2845 "Negative local idx");
2846 }
2847 }
2848
2849 } else {
2850
2851 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
2852 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
2853 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
2854 }
2855
2856 if (verb > QUIET) {
2858 "WORLD", Sev::inform,
2859 "Removed DOFs from problem %s dofs [%d / %d (before %d / %d) global]",
2860 prb_ptr->getName().c_str(), prb_ptr->getNbDofsRow(),
2861 prb_ptr->getNbDofsCol(), nb_init_row_dofs, nb_init_col_dofs);
2862 MOFEM_LOG_C("SYNC", Sev::verbose,
2863 "Removed DOFs from problem %s dofs [ %d / %d "
2864 "(before %d / %d) local, %d / %d (before %d / %d)]",
2865 prb_ptr->getName().c_str(), prb_ptr->getNbLocalDofsRow(),
2866 prb_ptr->getNbLocalDofsCol(), nb_init_loc_row_dofs,
2867 nb_init_loc_row_dofs, prb_ptr->getNbGhostDofsRow(),
2868 prb_ptr->getNbGhostDofsCol(), nb_init_ghost_row_dofs,
2869 nb_init_ghost_col_dofs);
2870 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
2871 }
2872
2874}
if(!static_cast< bool >(ifstream(param_file)))
#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 > numeredRowDofsPtr
store DOFs on rows for this problem