v0.13.2
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. More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblem (const std::string name, const bool square_matrix, int verb=VERBOSE)
 build problem data structures More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblem (Problem *problem_ptr, const bool square_matrix, int verb=VERBOSE)
 build problem data structures More...
 
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) More...
 
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) More...
 
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 More...
 
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 More...
 
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 More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionSimpleProblem (const std::string name, int verb=VERBOSE)
 partition problem dofs More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionProblem (const std::string name, int verb=VERBOSE)
 partition problem dofs (collective) More...
 
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 More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofs (const std::string name, int verb=VERBOSE)
 determine ghost nodes More...
 
MoFEMErrorCode MoFEM::ProblemsManager::partitionGhostDofsOnDistributedMesh (const std::string name, int verb=VERBOSE)
 determine ghost nodes on distributed meshes More...
 
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 More...
 
MoFEMErrorCode MoFEM::ProblemsManager::getProblemElementsLayout (const std::string name, const std::string &fe_name, PetscLayout *layout) const
 Get layout of elements in the problem. More...
 
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. More...
 
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. More...
 

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. More...
 

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<IS> *add_prb_is[] = {&cmp_prb_data->rowIs, &cmp_prb_data->colIs};
1288
1289 // Get local indices counter
1290 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1291 &out_problem_it->nbLocDofsCol};
1292 // Get global indices counter
1293 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1294
1295 // Set number of ghost nodes to zero
1296 {
1297 out_problem_it->nbDofsRow = 0;
1298 out_problem_it->nbDofsCol = 0;
1299 out_problem_it->nbLocDofsRow = 0;
1300 out_problem_it->nbLocDofsCol = 0;
1301 out_problem_it->nbGhostDofsRow = 0;
1302 out_problem_it->nbGhostDofsCol = 0;
1303 }
1304 int nb_dofs_reserve[] = {0, 0};
1305
1306 // Loop over rows and columns in the main problem and sub-problems
1307 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1308 // cerr << "SS " << ss << endl;
1309 // cerr << add_prb[ss]->size() << endl;
1310 // cerr << add_prb_ptr[ss]->size() << endl;
1311 add_prb_ptr[ss]->reserve(add_prb[ss]->size());
1312 add_prb_is[ss]->reserve(add_prb[ss]->size());
1313 for (std::vector<std::string>::const_iterator vit = add_prb[ss]->begin();
1314 vit != add_prb[ss]->end(); vit++) {
1315 ProblemByName::iterator prb_it = problems_by_name.find(*vit);
1316 if (prb_it == problems_by_name.end()) {
1317 SETERRQ1(
1318 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1319 "problem with name < %s > not defined (top tip check spelling)",
1320 vit->c_str());
1321 }
1322 add_prb_ptr[ss]->push_back(&*prb_it);
1323 // set number of dofs on rows and columns
1324 if (ss == 0) {
1325 // row
1326 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsRow();
1327 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsRow();
1328 nb_dofs_reserve[ss] +=
1329 add_prb_ptr[ss]->back()->numeredRowDofsPtr->size();
1330 } else {
1331 // column
1332 *nb_dofs[ss] += add_prb_ptr[ss]->back()->getNbDofsCol();
1333 *nb_local_dofs[ss] += add_prb_ptr[ss]->back()->getNbLocalDofsCol();
1334 nb_dofs_reserve[ss] +=
1335 add_prb_ptr[ss]->back()->numeredColDofsPtr->size();
1336 }
1337 }
1338 }
1339 // if squre problem, rows and columns are the same
1340 if (square_matrix) {
1341 add_prb_ptr[1]->reserve(add_prb_ptr[0]->size());
1342 add_prb_is[1]->reserve(add_prb_ptr[0]->size());
1343 out_problem_it->numeredColDofsPtr = out_problem_it->numeredRowDofsPtr;
1344 *nb_dofs[1] = *nb_dofs[0];
1345 *nb_local_dofs[1] = *nb_local_dofs[0];
1346 }
1347
1348 // reserve memory for dofs
1349 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array[2];
1350 // Reserve memory
1351 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1352 dofs_array[ss] = boost::make_shared<std::vector<NumeredDofEntity>>();
1353 dofs_array[ss]->reserve(nb_dofs_reserve[ss]);
1354 if (!ss)
1355 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array[ss]);
1356 else
1357 out_problem_it->getColDofsSequence()->emplace_back(dofs_array[ss]);
1358 }
1359
1360 // Push back DOFs
1361 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1362 NumeredDofEntity_multiIndex::index<PetscGlobalIdx_mi_tag>::type::iterator
1363 dit,
1364 hi_dit;
1365 int shift_glob = 0;
1366 int shift_loc = 0;
1367 for (unsigned int pp = 0; pp != add_prb_ptr[ss]->size(); pp++) {
1368 PetscInt *dofs_out_idx_ptr;
1369 int nb_local_dofs = (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1370 CHKERR PetscMalloc(nb_local_dofs * sizeof(int), &dofs_out_idx_ptr);
1371 if (ss == 0) {
1372 dit = (*add_prb_ptr[ss])[pp]
1373 ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1374 .begin();
1375 hi_dit = (*add_prb_ptr[ss])[pp]
1376 ->numeredRowDofsPtr->get<PetscGlobalIdx_mi_tag>()
1377 .end();
1378 } else {
1379 dit = (*add_prb_ptr[ss])[pp]
1380 ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1381 .begin();
1382 hi_dit = (*add_prb_ptr[ss])[pp]
1383 ->numeredColDofsPtr->get<PetscGlobalIdx_mi_tag>()
1384 .end();
1385 }
1386 int is_nb = 0;
1387 for (; dit != hi_dit; dit++) {
1388 const BitRefLevel &prb_bit = out_problem_it->getBitRefLevel();
1389 const BitRefLevel &prb_mask = out_problem_it->getBitRefLevelMask();
1390 const BitRefLevel &dof_bit = dit->get()->getBitRefLevel();
1391 if ((dof_bit & prb_bit).none() || ((dof_bit & prb_mask) != dof_bit))
1392 continue;
1393 const int rank = m_field.get_comm_rank();
1394 const int part = dit->get()->getPart();
1395 const int glob_idx = shift_glob + dit->get()->getPetscGlobalDofIdx();
1396 const int loc_idx =
1397 (part == rank) ? (shift_loc + dit->get()->getPetscLocalDofIdx())
1398 : -1;
1399 dofs_array[ss]->emplace_back(dit->get()->getDofEntityPtr(), glob_idx,
1400 glob_idx, loc_idx, part);
1401 if (part == rank) {
1402 dofs_out_idx_ptr[is_nb++] = glob_idx;
1403 }
1404 }
1405 if (is_nb > nb_local_dofs) {
1406 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1407 "Data inconsistency");
1408 }
1409 IS is;
1410 CHKERR ISCreateGeneral(m_field.get_comm(), is_nb, dofs_out_idx_ptr,
1411 PETSC_OWN_POINTER, &is);
1412 (*add_prb_is[ss]).push_back(is);
1413 if (ss == 0) {
1414 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsRow();
1415 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsRow();
1416 } else {
1417 shift_glob += (*add_prb_ptr[ss])[pp]->getNbDofsCol();
1418 shift_loc += (*add_prb_ptr[ss])[pp]->getNbLocalDofsCol();
1419 }
1420 if (square_matrix) {
1421 (*add_prb_ptr[1]).push_back((*add_prb_ptr[0])[pp]);
1422 (*add_prb_is[1]).push_back(is);
1423 CHKERR PetscObjectReference((PetscObject)is);
1424 }
1425 }
1426 }
1427
1428 if ((*add_prb_is[1]).size() != (*add_prb_is[0]).size()) {
1429 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1430 }
1431
1432 // Insert DOFs to problem multi-index
1433 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1434 auto hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->end()
1435 : out_problem_it->numeredColDofsPtr->end();
1436 for (auto &v : *dofs_array[ss])
1437 hint = (ss == 0) ? out_problem_it->numeredRowDofsPtr->emplace_hint(
1438 hint, dofs_array[ss], &v)
1439 : out_problem_it->numeredColDofsPtr->emplace_hint(
1440 hint, dofs_array[ss], &v);
1441 }
1442
1443 // Compress DOFs
1444 *nb_dofs[0] = 0;
1445 *nb_dofs[1] = 0;
1446 *nb_local_dofs[0] = 0;
1447 *nb_local_dofs[1] = 0;
1448 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ss++) {
1449
1450 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr;
1451 if (ss == 0) {
1452 dofs_ptr = out_problem_it->numeredRowDofsPtr;
1453 } else {
1454 dofs_ptr = out_problem_it->numeredColDofsPtr;
1455 }
1456 NumeredDofEntityByUId::iterator dit, hi_dit;
1457 dit = dofs_ptr->get<Unique_mi_tag>().begin();
1458 hi_dit = dofs_ptr->get<Unique_mi_tag>().end();
1459 std::vector<int> idx;
1460 idx.reserve(std::distance(dit, hi_dit));
1461 // set dofs in order entity and dof number on entity
1462 for (; dit != hi_dit; dit++) {
1463 if (dit->get()->getPart() == (unsigned int)m_field.get_comm_rank()) {
1464 bool success = dofs_ptr->get<Unique_mi_tag>().modify(
1465 dit, NumeredDofEntity_local_idx_change((*nb_local_dofs[ss])++));
1466 if (!success) {
1467 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1468 "modification unsuccessful");
1469 }
1470 idx.push_back(dit->get()->getPetscGlobalDofIdx());
1471 } else {
1472 if (dit->get()->getPetscLocalDofIdx() != -1) {
1473 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1474 "local index should be negative");
1475 }
1476 }
1477 }
1478 if (square_matrix) {
1479 *nb_local_dofs[1] = *nb_local_dofs[0];
1480 }
1481
1482 // set new dofs mapping
1483 IS is;
1484 CHKERR ISCreateGeneral(m_field.get_comm(), idx.size(), &*idx.begin(),
1485 PETSC_USE_POINTER, &is);
1486 CHKERR ISGetSize(is, nb_dofs[ss]);
1487 if (square_matrix) {
1488 *nb_dofs[1] = *nb_dofs[0];
1489 }
1490
1491 AO ao;
1492 CHKERR AOCreateMappingIS(is, PETSC_NULL, &ao);
1493 for (unsigned int pp = 0; pp != (*add_prb_is[ss]).size(); pp++)
1494 CHKERR AOApplicationToPetscIS(ao, (*add_prb_is[ss])[pp]);
1495
1496 // Set DOFs numeration
1497 {
1498 std::vector<int> idx_new;
1499 idx_new.reserve(dofs_ptr->size());
1500 for (NumeredDofEntityByUId::iterator dit =
1501 dofs_ptr->get<Unique_mi_tag>().begin();
1502 dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1503 idx_new.push_back(dit->get()->getPetscGlobalDofIdx());
1504 }
1505 // set new global dofs numeration
1506 IS is_new;
1507 CHKERR ISCreateGeneral(m_field.get_comm(), idx_new.size(),
1508 &*idx_new.begin(), PETSC_USE_POINTER, &is_new);
1509 CHKERR AOApplicationToPetscIS(ao, is_new);
1510 // set global indices to multi-index
1511 std::vector<int>::iterator vit = idx_new.begin();
1512 for (NumeredDofEntityByUId::iterator dit =
1513 dofs_ptr->get<Unique_mi_tag>().begin();
1514 dit != dofs_ptr->get<Unique_mi_tag>().end(); dit++) {
1515 bool success =
1516 dofs_ptr->modify(dit, NumeredDofEntity_part_and_glob_idx_change(
1517 dit->get()->getPart(), *(vit++)));
1518 if (!success) {
1519 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1520 "modification unsuccessful");
1521 }
1522 }
1523 CHKERR ISDestroy(&is_new);
1524 }
1525 CHKERR ISDestroy(&is);
1526 CHKERR AODestroy(&ao);
1527 }
1528
1529 CHKERR printPartitionedProblem(&*out_problem_it, verb);
1530 CHKERR debugPartitionedProblem(&*out_problem_it, verb);
1531
1532 // Inidcate that porble has been build
1535
1537}
#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()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
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.
Definition: LogManager.hpp:352
@ NOISY
Definition: definitions.h:211
@ VERBOSE
Definition: definitions.h:209
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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.
Definition: LogManager.hpp:308
uint128_t UId
Unique Id.
Definition: Types.hpp:31
const double r
rate factor
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 fe_ent_ptr = m_field.get_ents_finite_elements();
319 auto ents_field_ptr = m_field.get_field_ents();
320 auto dofs_field_ptr = m_field.get_dofs();
322 PetscLogEventBegin(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
323
324 // clear data structures
325 CHKERR m_field.clear_problem(problem_ptr->getName());
326
328
329 if (problem_ptr->getBitRefLevel().none())
330 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
331 "problem <%s> refinement level not set",
332 problem_ptr->getName().c_str());
333
334 int loop_size = 2;
335 if (square_matrix) {
336 loop_size = 1;
337 problem_ptr->numeredColDofsPtr = problem_ptr->numeredRowDofsPtr;
338 } else if (problem_ptr->numeredColDofsPtr == problem_ptr->numeredRowDofsPtr) {
339 problem_ptr->numeredColDofsPtr =
340 boost::shared_ptr<NumeredDofEntity_multiIndex>(
342 }
343
344 const BitRefLevel &prb_bit = problem_ptr->getBitRefLevel();
345 const BitRefLevel &prb_mask = problem_ptr->getBitRefLevelMask();
346
347 // // get rows and cols dofs view based on data on elements
348 DofEntity_multiIndex_active_view dofs_rows, dofs_cols;
349
350 // Add DOFs to problem by visiting all elements and adding DOFs from
351 // elements to the problem
352 if (buildProblemFromFields == PETSC_FALSE) {
353
354 auto make_rows_and_cols_view = [&](auto &dofs_rows, auto &dofs_cols) {
355 auto ents_field_ptr = m_field.get_field_ents();
356 auto adjacencies_ptr = m_field.get_ents_elements_adjacency();
358 for (auto it = ents_field_ptr->begin(); it != ents_field_ptr->end();
359 ++it) {
360
361 const auto uid = (*it)->getLocalUniqueId();
362
363 auto r = adjacencies_ptr->get<Unique_mi_tag>().equal_range(uid);
364 for (auto lo = r.first; lo != r.second; ++lo) {
365
366 if ((lo->getBitFEId() & problem_ptr->getBitFEId()).any()) {
367 std::array<bool, 2> row_col = {false, false};
368
369 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
370
371 // if entity is not problem refinement level
372 if ((fe_bit & prb_bit).any() && (fe_bit & prb_mask) == fe_bit) {
373
374 auto add_to_view = [&](auto dofs, auto &view, auto rc) {
375 auto dit = dofs->lower_bound(uid);
376 decltype(dit) hi_dit;
377 if (dit != dofs->end()) {
378 hi_dit = dofs->upper_bound(
379 uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
380 view.insert(dit, hi_dit);
381 row_col[rc] = true;
382 }
383 };
384
385 if ((lo->entFePtr->getBitFieldIdRow() & (*it)->getId()).any())
386 add_to_view(dofs_field_ptr, dofs_rows, ROW);
387
388 if (!square_matrix)
389 if ((lo->entFePtr->getBitFieldIdCol() & (*it)->getId()).any())
390 add_to_view(dofs_field_ptr, dofs_cols, COL);
391
392 if (square_matrix && row_col[ROW])
393 break;
394 else if (row_col[ROW] && row_col[COL])
395 break;
396 }
397 }
398 }
399 }
401 };
402
403 CHKERR make_rows_and_cols_view(dofs_rows, dofs_cols);
404 }
405
406 // Add DOFS to the problem by searching all the filedes, and adding to
407 // problem owned or shared DOFs
408 if (buildProblemFromFields == PETSC_TRUE) {
409 // Get fields IDs on elements
410 BitFieldId fields_ids_row, fields_ids_col;
411 for (auto fit = fe_ptr->begin(); fit != fe_ptr->end(); fit++) {
412 if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
413 fields_ids_row |= fit->get()->getBitFieldIdRow();
414 fields_ids_col |= fit->get()->getBitFieldIdCol();
415 }
416 }
417 // Get fields DOFs
418 for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
419 if ((fit->get()->getId() & (fields_ids_row | fields_ids_col)).any()) {
420
421 auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(
422 FieldEntity::getLoBitNumberUId((*fit)->getBitNumber()));
423 auto hi_dit = dofs_field_ptr->get<Unique_mi_tag>().upper_bound(
424 FieldEntity::getHiBitNumberUId((*fit)->getBitNumber()));
425
426 for (; dit != hi_dit; dit++) {
427
428 const int owner_proc = dit->get()->getOwnerProc();
429 if (owner_proc != m_field.get_comm_rank()) {
430 const unsigned char pstatus = dit->get()->getPStatus();
431 if (pstatus == 0) {
432 continue;
433 }
434 }
435
436 const auto &dof_bit = (*dit)->getBitRefLevel();
437 // if entity is not problem refinement level
438 if ((dof_bit & prb_bit).any() && (dof_bit & prb_mask) == dof_bit) {
439
440 if ((fit->get()->getId() & fields_ids_row).any()) {
441 dofs_rows.insert(*dit);
442 }
443 if (!square_matrix) {
444 if ((fit->get()->getId() & fields_ids_col).any()) {
445 dofs_cols.insert(*dit);
446 }
447 }
448
449 }
450 }
451 }
452 }
453 }
454
456 // Get fields IDs on elements
457 BitFieldId fields_ids_row, fields_ids_col;
458 BitFieldId *fields_ids[2] = {&fields_ids_row, &fields_ids_col};
459 for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
460 fit != fe_ptr->end(); fit++) {
461 if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
462 fields_ids_row |= fit->get()->getBitFieldIdRow();
463 fields_ids_col |= fit->get()->getBitFieldIdCol();
464 }
465 }
466
467 DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
468 for (int ss = 0; ss != ((square_matrix) ? 1 : 2); ++ss) {
469 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
470 hi_miit;
471 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
472 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
473 Range ents_to_synchronise;
474 for (; miit != hi_miit; ++miit) {
475 if (miit->get()->getEntDofIdx() != 0)
476 continue;
477 ents_to_synchronise.insert(miit->get()->getEnt());
478 }
479 Range tmp_ents = ents_to_synchronise;
480 CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
481 ents_to_synchronise, verb);
482 ents_to_synchronise = subtract(ents_to_synchronise, tmp_ents);
483 for (auto fit = fields_ptr->begin(); fit != fields_ptr->end(); fit++) {
484 if ((fit->get()->getId() & *fields_ids[ss]).any()) {
485 const auto bit_number = (*fit)->getBitNumber();
486 for (Range::pair_iterator pit = ents_to_synchronise.pair_begin();
487 pit != ents_to_synchronise.pair_end(); ++pit) {
488 const auto f = pit->first;
489 const auto s = pit->second;
490 const auto lo_uid =
492 const auto hi_uid =
494
495 auto dit = dofs_field_ptr->get<Unique_mi_tag>().lower_bound(lo_uid);
496 auto hi_dit =
497 dofs_field_ptr->get<Unique_mi_tag>().upper_bound(hi_uid);
498
499 dofs_ptr[ss]->insert(dit, hi_dit);
500 }
501 }
502 }
503 }
504 }
505
506 // add dofs for rows and cols and set ownership
507 DofEntity_multiIndex_active_view *dofs_ptr[] = {&dofs_rows, &dofs_cols};
508 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_ptr[] = {
509 problem_ptr->numeredRowDofsPtr, problem_ptr->numeredColDofsPtr};
510 int *nbdof_ptr[] = {&problem_ptr->nbDofsRow, &problem_ptr->nbDofsCol};
511 int *local_nbdof_ptr[] = {&problem_ptr->nbLocDofsRow,
512 &problem_ptr->nbLocDofsCol};
513 int *ghost_nbdof_ptr[] = {&problem_ptr->nbGhostDofsRow,
514 &problem_ptr->nbGhostDofsCol};
515 for (int ss = 0; ss < 2; ss++) {
516 *(nbdof_ptr[ss]) = 0;
517 *(local_nbdof_ptr[ss]) = 0;
518 *(ghost_nbdof_ptr[ss]) = 0;
519 }
520
521 // Loop over dofs on rows and columns and add to multi-indices in dofs
522 // problem structure, set partition for each dof
523 int nb_local_dofs[] = {0, 0};
524 for (int ss = 0; ss < loop_size; ss++) {
525 DofEntity_multiIndex_active_view::nth_index<1>::type::iterator miit,
526 hi_miit;
527 miit = dofs_ptr[ss]->get<1>().lower_bound(1);
528 hi_miit = dofs_ptr[ss]->get<1>().upper_bound(1);
529 for (; miit != hi_miit; miit++) {
530 int owner_proc = (*miit)->getOwnerProc();
531 if (owner_proc == m_field.get_comm_rank()) {
532 nb_local_dofs[ss]++;
533 }
534 }
535 }
536
537 // get layout
538 int start_ranges[2], end_ranges[2];
539 for (int ss = 0; ss != loop_size; ss++) {
540 PetscLayout layout;
541 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
542 CHKERR PetscLayoutSetBlockSize(layout, 1);
543 CHKERR PetscLayoutSetLocalSize(layout, nb_local_dofs[ss]);
544 CHKERR PetscLayoutSetUp(layout);
545 CHKERR PetscLayoutGetSize(layout, &*nbdof_ptr[ss]);
546 CHKERR PetscLayoutGetRange(layout, &start_ranges[ss], &end_ranges[ss]);
547 CHKERR PetscLayoutDestroy(&layout);
548 }
549 if (square_matrix) {
550 nbdof_ptr[1] = nbdof_ptr[0];
551 nb_local_dofs[1] = nb_local_dofs[0];
552 start_ranges[1] = start_ranges[0];
553 end_ranges[1] = end_ranges[0];
554 }
555
556 // if(sizeof(UId) != SIZEOFUID) {
557 // SETERRQ2(
558 // PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,
559 // "check size of UId, size of UId is %u != %u",
560 // sizeof(UId),SIZEOFUID
561 // );
562 // }
563
564 // set local and global indices on own dofs
565 const size_t idx_data_type_size = sizeof(IdxDataType);
566 const size_t data_block_size = idx_data_type_size / sizeof(int);
567
568 if (sizeof(IdxDataType) % sizeof(int)) {
569 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
570 }
571 std::vector<std::vector<IdxDataType>> ids_data_packed_rows(
572 m_field.get_comm_size()),
573 ids_data_packed_cols(m_field.get_comm_size());
574
575 // Loop over dofs on this processor and prepare those dofs to send on
576 // another proc
577 for (int ss = 0; ss != loop_size; ++ss) {
578
579 DofEntity_multiIndex_active_view::nth_index<0>::type::iterator miit,
580 hi_miit;
581 hi_miit = dofs_ptr[ss]->get<0>().end();
582
583 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
584 boost::shared_ptr<std::vector<NumeredDofEntity>>(
585 new std::vector<NumeredDofEntity>());
586 int nb_dofs_to_add = 0;
587 miit = dofs_ptr[ss]->get<0>().begin();
588 for (; miit != hi_miit; ++miit) {
589 // Only set global idx for dofs on this processor part
590 if (!(miit->get()->getActive()))
591 continue;
592 ++nb_dofs_to_add;
593 }
594 dofs_array->reserve(nb_dofs_to_add);
595 if (ss == 0) {
596 problem_ptr->getRowDofsSequence()->push_back(dofs_array);
597 } else {
598 problem_ptr->getColDofsSequence()->push_back(dofs_array);
599 }
600
601 int &local_idx = *local_nbdof_ptr[ss];
602 miit = dofs_ptr[ss]->get<0>().begin();
603 for (; miit != hi_miit; ++miit) {
604
605 // Only set global idx for dofs on this processor part
606 if (!(miit->get()->getActive()))
607 continue;
608
609 dofs_array->emplace_back(*miit);
610
611 int owner_proc = dofs_array->back().getOwnerProc();
612 if (owner_proc < 0) {
613 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
614 "data inconsistency");
615 }
616
617 if (owner_proc != m_field.get_comm_rank()) {
618 dofs_array->back().pArt = owner_proc;
619 dofs_array->back().dofIdx = -1;
620 dofs_array->back().petscGloablDofIdx = -1;
621 dofs_array->back().petscLocalDofIdx = -1;
622 } else {
623
624 // Set part and indexes
625 int glob_idx = start_ranges[ss] + local_idx;
626 dofs_array->back().pArt = owner_proc;
627 dofs_array->back().dofIdx = glob_idx;
628 dofs_array->back().petscGloablDofIdx = glob_idx;
629 dofs_array->back().petscLocalDofIdx = local_idx;
630 local_idx++;
631
632 unsigned char pstatus = dofs_array->back().getPStatus();
633 // check id dof is shared, if that is a case global idx added to data
634 // structure and is sended to other processors
635 if (pstatus > 0) {
636
637 for (int proc = 0;
638 proc < MAX_SHARING_PROCS &&
639 -1 != dofs_array->back().getSharingProcsPtr()[proc];
640 proc++) {
641 // make it different for rows and columns when needed
642 if (ss == 0) {
643 ids_data_packed_rows[dofs_array->back()
644 .getSharingProcsPtr()[proc]]
645 .emplace_back(dofs_array->back().getGlobalUniqueId(),
646 glob_idx);
647 } else {
648 ids_data_packed_cols[dofs_array->back()
649 .getSharingProcsPtr()[proc]]
650 .emplace_back(dofs_array->back().getGlobalUniqueId(),
651 glob_idx);
652 }
653 if (!(pstatus & PSTATUS_MULTISHARED)) {
654 break;
655 }
656 }
657 }
658 }
659 }
660
661 auto hint = numered_dofs_ptr[ss]->end();
662 for (auto &v : *dofs_array)
663 hint = numered_dofs_ptr[ss]->emplace_hint(hint, dofs_array, &v);
664 }
665 if (square_matrix) {
666 local_nbdof_ptr[1] = local_nbdof_ptr[0];
667 }
668
669 int nsends_rows = 0, nsends_cols = 0;
670 // Non zero lengths[i] represent a message to i of length lengths[i].
671 std::vector<int> lengths_rows(m_field.get_comm_size()),
672 lengths_cols(m_field.get_comm_size());
673 lengths_rows.clear();
674 lengths_cols.clear();
675 for (int proc = 0; proc < m_field.get_comm_size(); proc++) {
676 lengths_rows[proc] = ids_data_packed_rows[proc].size() * data_block_size;
677 lengths_cols[proc] = ids_data_packed_cols[proc].size() * data_block_size;
678 if (!ids_data_packed_rows[proc].empty())
679 nsends_rows++;
680 if (!ids_data_packed_cols[proc].empty())
681 nsends_cols++;
682 }
683
684 MPI_Status *status;
685 CHKERR PetscMalloc1(m_field.get_comm_size(), &status);
686
687 // Do rows
688 int nrecvs_rows; // number of messages received
689 int *onodes_rows; // list of node-ids from which messages are expected
690 int *olengths_rows; // corresponding message lengths
691 int **rbuf_row; // must bee freed by user
692
693 // make sure it is a PETSc comm
694 MPI_Comm comm;
695 CHKERR PetscCommDuplicate(m_field.get_comm(), &comm, NULL);
696
697 {
698
699 // rows
700
701 // Computes the number of messages a node expects to receive
702 CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_rows[0],
703 &nrecvs_rows);
704 // std::cerr << nrecvs_rows << std::endl;
705
706 // Computes info about messages that a MPI-node will receive, including
707 // (from-id,length) pairs for each message.
708 CHKERR PetscGatherMessageLengths(comm, nsends_rows, nrecvs_rows,
709 &lengths_rows[0], &onodes_rows,
710 &olengths_rows);
711
712 // Gets a unique new tag from a PETSc communicator. All processors that
713 // share the communicator MUST call this routine EXACTLY the same number
714 // of times. This tag should only be used with the current objects
715 // communicator; do NOT use it with any other MPI communicator.
716 int tag_row;
717 CHKERR PetscCommGetNewTag(comm, &tag_row);
718
719 // Allocate a buffer sufficient to hold messages of size specified in
720 // olengths. And post Irecvs on these buffers using node info from onodes
721 MPI_Request *r_waits_row; // must bee freed by user
722 // rbuf has a pointers to messeges. It has size of of nrecvs (number of
723 // messages) +1. In the first index a block is allocated,
724 // such that rbuf[i] = rbuf[i-1]+olengths[i-1].
725
726 CHKERR PetscPostIrecvInt(comm, tag_row, nrecvs_rows, onodes_rows,
727 olengths_rows, &rbuf_row, &r_waits_row);
728 CHKERR PetscFree(onodes_rows);
729
730 MPI_Request *s_waits_row; // status of sens messages
731 CHKERR PetscMalloc1(nsends_rows, &s_waits_row);
732
733 // Send messeges
734 for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
735 if (!lengths_rows[proc])
736 continue; // no message to send to this proc
737 CHKERR MPI_Isend(&(ids_data_packed_rows[proc])[0], // buffer to send
738 lengths_rows[proc], // message length
739 MPIU_INT, proc, // to proc
740 tag_row, comm, s_waits_row + kk);
741 kk++;
742 }
743
744 if (nrecvs_rows) {
745 CHKERR MPI_Waitall(nrecvs_rows, r_waits_row, status);
746 }
747 if (nsends_rows) {
748 CHKERR MPI_Waitall(nsends_rows, s_waits_row, status);
749 }
750
751 CHKERR PetscFree(r_waits_row);
752 CHKERR PetscFree(s_waits_row);
753 }
754
755 // cols
756 int nrecvs_cols = nrecvs_rows;
757 int *olengths_cols = olengths_rows;
758 PetscInt **rbuf_col = rbuf_row;
759 if (!square_matrix) {
760
761 // Computes the number of messages a node expects to receive
762 CHKERR PetscGatherNumberOfMessages(comm, NULL, &lengths_cols[0],
763 &nrecvs_cols);
764
765 // Computes info about messages that a MPI-node will receive, including
766 // (from-id,length) pairs for each message.
767 int *onodes_cols;
768 CHKERR PetscGatherMessageLengths(comm, nsends_cols, nrecvs_cols,
769 &lengths_cols[0], &onodes_cols,
770 &olengths_cols);
771
772 // Gets a unique new tag from a PETSc communicator.
773 int tag_col;
774 CHKERR PetscCommGetNewTag(comm, &tag_col);
775
776 // Allocate a buffer sufficient to hold messages of size specified in
777 // olengths. And post Irecvs on these buffers using node info from onodes
778 MPI_Request *r_waits_col; // must bee freed by user
779 CHKERR PetscPostIrecvInt(comm, tag_col, nrecvs_cols, onodes_cols,
780 olengths_cols, &rbuf_col, &r_waits_col);
781 CHKERR PetscFree(onodes_cols);
782
783 MPI_Request *s_waits_col; // status of sens messages
784 CHKERR PetscMalloc1(nsends_cols, &s_waits_col);
785
786 // Send messages
787 for (int proc = 0, kk = 0; proc < m_field.get_comm_size(); proc++) {
788 if (!lengths_cols[proc])
789 continue; // no message to send to this proc
790 CHKERR MPI_Isend(&(ids_data_packed_cols[proc])[0], // buffer to send
791 lengths_cols[proc], // message length
792 MPIU_INT, proc, // to proc
793 tag_col, comm, s_waits_col + kk);
794 kk++;
795 }
796
797 if (nrecvs_cols) {
798 CHKERR MPI_Waitall(nrecvs_cols, r_waits_col, status);
799 }
800 if (nsends_cols) {
801 CHKERR MPI_Waitall(nsends_cols, s_waits_col, status);
802 }
803
804 CHKERR PetscFree(r_waits_col);
805 CHKERR PetscFree(s_waits_col);
806 }
807
808 CHKERR PetscCommDestroy(&comm);
809 CHKERR PetscFree(status);
810
811 DofEntity_multiIndex_global_uid_view dofs_glob_uid_view;
812 auto hint = dofs_glob_uid_view.begin();
813 for (auto dof : *m_field.get_dofs())
814 dofs_glob_uid_view.emplace_hint(hint, dof);
815
816 // set values received from other processors
817 for (int ss = 0; ss != loop_size; ++ss) {
818
819 int nrecvs;
820 int *olengths;
821 int **data_procs;
822 if (ss == 0) {
823 nrecvs = nrecvs_rows;
824 olengths = olengths_rows;
825 data_procs = rbuf_row;
826 } else {
827 nrecvs = nrecvs_cols;
828 olengths = olengths_cols;
829 data_procs = rbuf_col;
830 }
831
832 UId uid;
833 for (int kk = 0; kk != nrecvs; ++kk) {
834 int len = olengths[kk];
835 int *data_from_proc = data_procs[kk];
836 for (int dd = 0; dd < len; dd += data_block_size) {
837 uid = IdxDataTypePtr(&data_from_proc[dd]).getUId();
838 auto ddit = dofs_glob_uid_view.find(uid);
839
840 if (PetscUnlikely(ddit == dofs_glob_uid_view.end())) {
841
842#ifndef NDEBUG
843 MOFEM_LOG("SELF", Sev::error)
844 << "DOF is shared or multishared between processors. For example "
845 "if order of field on given entity is set in inconsistently, "
846 "has different value on two processor, error such as this is "
847 "triggered";
848
849 MOFEM_LOG("SELF", Sev::error) << "UId " << uid << " is not found";
850 const auto owner_proc = FieldEntity::getOwnerFromUniqueId(uid);
851 MOFEM_LOG("SELF", Sev::error)
852 << "Problematic UId owner proc is " << owner_proc;
853 const auto uid_handle = FieldEntity::getHandleFromUniqueId(uid);
854 MOFEM_LOG("SELF", Sev::error)
855 << "Problematic UId entity owning processor handle is "
856 << uid_handle << " entity type "
857 << moab::CN::EntityTypeName(type_from_handle(uid_handle));
858 const auto uid_bit_number =
860 MOFEM_LOG("SELF", Sev::error)
861 << "Problematic UId field is "
862 << m_field.get_field_name(
863 field_bit_from_bit_number(uid_bit_number));
864
866 field_view.insert(ents_field_ptr->begin(), ents_field_ptr->end());
867 auto fe_it = field_view.find(FieldEntity::getGlobalUniqueIdCalculate(
868 owner_proc, uid_bit_number, uid_handle));
869 if (fe_it == field_view.end()) {
870 MOFEM_LOG("SELF", Sev::error)
871 << "Also, no field entity in database for given global UId";
872 } else {
873 MOFEM_LOG("SELF", Sev::error) << "Field entity in databse exist "
874 "(but have no DOF wih give UId";
875 MOFEM_LOG("SELF", Sev::error) << **fe_it;
876
877 // Save file with missing entity
878 auto error_file_name =
879 "error_with_missing_entity_" +
880 boost::lexical_cast<std::string>(m_field.get_comm_rank()) +
881 ".vtk";
882 MOFEM_LOG("SELF", Sev::error)
883 << "Look to file < " << error_file_name
884 << " > it contains entity with missing DOF.";
885
886 auto tmp_msh = get_temp_meshset_ptr(m_field.get_moab());
887 const auto local_fe_ent = (*fe_it)->getEnt();
888 CHKERR m_field.get_moab().add_entities(*tmp_msh, &local_fe_ent, 1);
889 CHKERR m_field.get_moab().write_file(error_file_name.c_str(), "VTK",
890 "", tmp_msh->get_ptr(), 1);
891 }
892#endif
893
894 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
895 "DOF with global UId not found (Compile code in Debug to "
896 "learn more about problem");
897 }
898
899 auto dit = numered_dofs_ptr[ss]->find((*ddit)->getLocalUniqueId());
900
901 if (dit != numered_dofs_ptr[ss]->end()) {
902
903 int global_idx = IdxDataTypePtr(&data_from_proc[dd]).getDofIdx();
904#ifndef NDEBUG
905 if (PetscUnlikely(global_idx < 0))
906 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
907 "received negative dof");
908#endif
909 bool success;
910 success = numered_dofs_ptr[ss]->modify(
911 dit, NumeredDofEntity_mofem_index_change(global_idx));
912
913#ifndef NDEBUG
914 if (PetscUnlikely(!success))
915 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
916 "modification unsuccessful");
917#endif
918 success = numered_dofs_ptr[ss]->modify(
919 dit, NumeredDofEntity_part_and_glob_idx_change((*dit)->getPart(),
920 global_idx));
921#ifndef NDEBUG
922 if (PetscUnlikely(!success))
923 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
924 "modification unsuccessful");
925#endif
926 };
927
928#ifndef NDEBUG
929 if (PetscUnlikely(ddit->get()->getPStatus() == 0)) {
930
931 // Dof is shared on this processor, however there is no element
932 // which have this dof. If DOF is not shared and received from other
933 // processor, but not marked as a shared on other that means that is
934 // data inconstancy and error should be thrown.
935
936 std::ostringstream zz;
937 zz << **ddit << std::endl;
938 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
939 "data inconsistency, dofs is not shared, but received "
940 "from other proc\n"
941 "%s",
942 zz.str().c_str());
943 }
944#endif
945 }
946 }
947 }
948
949 if (square_matrix) {
950 (problem_ptr->nbDofsCol) = (problem_ptr->nbDofsRow);
951 (problem_ptr->nbLocDofsCol) = (problem_ptr->nbLocDofsRow);
952 }
953
954 CHKERR PetscFree(olengths_rows);
955 CHKERR PetscFree(rbuf_row[0]);
956 CHKERR PetscFree(rbuf_row);
957 if (!square_matrix) {
958 CHKERR PetscFree(olengths_cols);
959 CHKERR PetscFree(rbuf_col[0]);
960 CHKERR PetscFree(rbuf_col);
961 }
962
963 if (square_matrix) {
964 if (numered_dofs_ptr[0]->size() != numered_dofs_ptr[1]->size()) {
965 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
966 "data inconsistency for square_matrix %d!=%d",
967 numered_dofs_ptr[0]->size(), numered_dofs_ptr[1]->size());
968 }
969 if (problem_ptr->numeredRowDofsPtr != problem_ptr->numeredColDofsPtr) {
970 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
971 "data inconsistency for square_matrix");
972 }
973 }
974
975 CHKERR printPartitionedProblem(problem_ptr, verb);
976 CHKERR debugPartitionedProblem(problem_ptr, verb);
977
978 PetscLogEventEnd(MOFEM_EVENT_ProblemsManager, 0, 0, 0, 0);
979
981}
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
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
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
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
virtual const EntFiniteElement_multiIndex * get_ents_finite_elements() const =0
Get the ents 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
auto f
Definition: HenckyOps.hpp:5
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:42
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1634
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
Definition: Templates.hpp:1650
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1626
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)
static auto getFieldBitNumberFromUniqueId(const UId uid)
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
remove_entities_from_problem.cpp.

Definition at line 983 of file ProblemsManager.cpp.

994 {
995 MoFEM::Interface &m_field = cOre;
996 auto problems_ptr = m_field.get_problems();
998
999 CHKERR m_field.clear_problem(out_name);
1000
1001 // get reference to all problems
1002 using ProblemByName = decltype(problems_ptr->get<Problem_mi_tag>());
1003 auto &problems_by_name =
1004 const_cast<ProblemByName &>(problems_ptr->get<Problem_mi_tag>());
1005
1006 // get iterators to out problem, i.e. build problem
1007 auto out_problem_it = problems_by_name.find(out_name);
1008 if (out_problem_it == problems_by_name.end()) {
1009 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1010 "subproblem with name < %s > not defined (top tip check spelling)",
1011 out_name.c_str());
1012 }
1013 // get iterator to main problem, i.e. out problem is subproblem of main
1014 // problem
1015 auto main_problem_it = problems_by_name.find(main_problem);
1016 if (main_problem_it == problems_by_name.end()) {
1017 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1018 "problem of subproblem with name < %s > not defined (top tip "
1019 "check spelling)",
1020 main_problem.c_str());
1021 }
1022
1023 // get dofs for row & columns for out problem,
1024 boost::shared_ptr<NumeredDofEntity_multiIndex> out_problem_dofs[] = {
1025 out_problem_it->numeredRowDofsPtr, out_problem_it->numeredColDofsPtr};
1026 // get dofs for row & columns for main problem
1027 boost::shared_ptr<NumeredDofEntity_multiIndex> main_problem_dofs[] = {
1028 main_problem_it->numeredRowDofsPtr, main_problem_it->numeredColDofsPtr};
1029 // get local indices counter
1030 int *nb_local_dofs[] = {&out_problem_it->nbLocDofsRow,
1031 &out_problem_it->nbLocDofsCol};
1032 // get global indices counter
1033 int *nb_dofs[] = {&out_problem_it->nbDofsRow, &out_problem_it->nbDofsCol};
1034
1035 // set number of ghost nodes to zero
1036 {
1037 out_problem_it->nbGhostDofsRow = 0;
1038 out_problem_it->nbGhostDofsCol = 0;
1039 }
1040
1041 // put rows & columns field names in array
1042 std::vector<std::string> fields[] = {fields_row, fields_col};
1043 const map<std::string, boost::shared_ptr<Range>> *entityMap[] = {
1044 entityMapRow, entityMapCol};
1045
1046 // make data structure fos sub-problem data
1047 out_problem_it->subProblemData =
1048 boost::make_shared<Problem::SubProblemData>();
1049
1050 // Loop over rows and columns
1051 for (int ss = 0; ss != (square_matrix ? 1 : 2); ++ss) {
1052
1053 // reset dofs and columns counters
1054 (*nb_local_dofs[ss]) = 0;
1055 (*nb_dofs[ss]) = 0;
1056 // clear arrays
1057 out_problem_dofs[ss]->clear();
1058
1059 // If DOFs are cleared clear finite elements too.
1060 out_problem_it->numeredFiniteElementsPtr->clear();
1061
1062 // get dofs by field name and insert them in out problem multi-indices
1063 for (auto field : fields[ss]) {
1064
1065 // Following reserve memory in sequences, only two allocations are here,
1066 // once for array of objects, next for array of shared pointers
1067
1068 // aliased sequence of pointer is killed with element
1069 boost::shared_ptr<std::vector<NumeredDofEntity>> dofs_array =
1070 boost::make_shared<std::vector<NumeredDofEntity>>();
1071 // reserve memory for field dofs
1072 if (!ss)
1073 out_problem_it->getRowDofsSequence()->emplace_back(dofs_array);
1074 else
1075 out_problem_it->getColDofsSequence()->emplace_back(dofs_array);
1076
1077 // create elements objects
1078 auto bit_number = m_field.get_field_bit_number(field);
1079
1080 auto add_dit_to_dofs_array = [&](auto &dit) {
1081 if (dit->get()->getPetscGlobalDofIdx() >= 0)
1082 dofs_array->emplace_back(
1083 dit->get()->getDofEntityPtr(), dit->get()->getPetscGlobalDofIdx(),
1084 dit->get()->getPetscGlobalDofIdx(),
1085 dit->get()->getPetscLocalDofIdx(), dit->get()->getPart());
1086 };
1087
1088 auto get_dafult_dof_range = [&]() {
1089 auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1090 FieldEntity::getLoBitNumberUId(bit_number));
1091 auto hi_dit = main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1092 FieldEntity::getHiBitNumberUId(bit_number));
1093 return std::make_pair(dit, hi_dit);
1094 };
1095
1096 if (entityMap[ss]) {
1097 auto mit = entityMap[ss]->find(field);
1098
1099 if (mit != entityMap[ss]->end()) {
1100 for (auto p = mit->second->pair_begin(); p != mit->second->pair_end();
1101 ++p) {
1102 const auto lo_ent = p->first;
1103 const auto hi_ent = p->second;
1104 auto dit = main_problem_dofs[ss]->get<Unique_mi_tag>().lower_bound(
1105 DofEntity::getLoFieldEntityUId(bit_number, lo_ent));
1106 auto hi_dit =
1107 main_problem_dofs[ss]->get<Unique_mi_tag>().upper_bound(
1108 DofEntity::getHiFieldEntityUId(bit_number, hi_ent));
1109 dofs_array->reserve(std::distance(dit, hi_dit));
1110 for (; dit != hi_dit; dit++) {
1111 add_dit_to_dofs_array(dit);
1112 }
1113 }
1114 } else {
1115 auto [dit, hi_dit] = get_dafult_dof_range();
1116 dofs_array->reserve(std::distance(dit, hi_dit));
1117 for (; dit != hi_dit; dit++)
1118 add_dit_to_dofs_array(dit);
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 2561 of file ProblemsManager.cpp.

2563 {
2564 MoFEM::Interface &m_field = cOre;
2565 const Problem *problem_ptr;
2566 const FiniteElement_multiIndex *fes_ptr;
2568
2569 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2570 CHKERR m_field.get_problem(prb_name, &problem_ptr);
2571 CHKERR m_field.get_finite_elements(&fes_ptr);
2572
2573 auto fe_miit = fes_ptr->get<FiniteElement_name_mi_tag>().find(fe_name);
2574 if (fe_miit != fes_ptr->get<FiniteElement_name_mi_tag>().end()) {
2575 auto fit =
2576 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
2578 0, (*fe_miit)->getFEUId()));
2579 auto hi_fe_it =
2580 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
2582 get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
2583 std::vector<EntityHandle> fe_vec;
2584 fe_vec.reserve(std::distance(fit, hi_fe_it));
2585 for (; fit != hi_fe_it; fit++)
2586 fe_vec.push_back(fit->get()->getEnt());
2587 CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2588 fe_vec.size());
2589 }
2590
2592}
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 2595 of file ProblemsManager.cpp.

2597 {
2598 MoFEM::Interface &m_field = cOre;
2599 const Problem *problem_ptr;
2601 CHKERR m_field.get_problem(name, &problem_ptr);
2602 CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2603 fe_name, layout);
2605}

◆ 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 1898 of file ProblemsManager.cpp.

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

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

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

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

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

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

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

◆ 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 3146 of file ProblemsManager.cpp.

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

Definition at line 2607 of file ProblemsManager.cpp.

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