v0.13.1
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, std::pair< EntityType, EntityType > > *entityMapRow=nullptr, const map< std::string, std::pair< EntityType, EntityType > > *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 1252 of file ProblemsManager.cpp.

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

103 {
104
105 MoFEM::Interface &m_field = cOre;
107 if (!(cOre.getBuildMoFEM() & (1 << 0)))
108 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
109 if (!(cOre.getBuildMoFEM() & (1 << 1)))
110 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
111 if (!(cOre.getBuildMoFEM() & (1 << 2)))
112 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
113 const Problem *problem_ptr;
114 CHKERR m_field.get_problem(name, &problem_ptr);
115 CHKERR buildProblem(const_cast<Problem *>(problem_ptr), square_matrix, verb);
116 cOre.getBuildMoFEM() |= 1 << 3; // It is assumed that user who uses this
117 // function knows what he is doing
119}
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 121 of file ProblemsManager.cpp.

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

305 {
306 MoFEM::Interface &m_field = cOre;
308
310 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
311 if (!((cOre.getBuildMoFEM()) & Core::BUILD_FE))
312 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
314 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
315
316 const Problem *problem_ptr;
317 CHKERR m_field.get_problem(name, &problem_ptr);
318 CHKERR buildProblemOnDistributedMesh(const_cast<Problem *>(problem_ptr),
319 square_matrix, verb);
320
323
325}
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 327 of file ProblemsManager.cpp.

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

Definition at line 997 of file ProblemsManager.cpp.

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

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

2565 {
2566 MoFEM::Interface &m_field = cOre;
2567 const Problem *problem_ptr;
2569
2570 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2571 CHKERR m_field.get_problem(prb_name, &problem_ptr);
2572 auto fit =
2573 problem_ptr->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
2574 .lower_bound(fe_name);
2575 auto hi_fe_it =
2576 problem_ptr->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
2577 .upper_bound(fe_name);
2578 std::vector<EntityHandle> fe_vec;
2579 fe_vec.reserve(std::distance(fit, hi_fe_it));
2580 for (; fit != hi_fe_it; fit++)
2581 fe_vec.push_back(fit->get()->getEnt());
2582 CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2583 fe_vec.size());
2585}
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements

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

2590 {
2591 MoFEM::Interface &m_field = cOre;
2592 const Problem *problem_ptr;
2594 CHKERR m_field.get_problem(name, &problem_ptr);
2595 CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2596 fe_name, layout);
2598}

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

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

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

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

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

Definition at line 2477 of file ProblemsManager.cpp.

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

533 {
534 MoFEM::Interface &m_field = cOre;
536
537 // get layout
538 int rstart, rend, nb_elems;
539 {
540 PetscLayout layout;
541 CHKERR PetscLayoutCreate(m_field.get_comm(), &layout);
542 CHKERR PetscLayoutSetBlockSize(layout, 1);
543 CHKERR PetscLayoutSetSize(layout, ents.size());
544 CHKERR PetscLayoutSetUp(layout);
545 CHKERR PetscLayoutGetSize(layout, &nb_elems);
546 CHKERR PetscLayoutGetRange(layout, &rstart, &rend);
547 CHKERR PetscLayoutDestroy(&layout);
548 if (verb >= VERBOSE) {
549 MOFEM_LOG("SYNC", Sev::inform)
550 << "Finite elements in problem: row lower " << rstart << " row upper "
551 << rend << " nb. elems " << nb_elems << " ( " << ents.size() << " )";
553 }
554 }
555
556 std::vector<EntityHandle> weight_ents;
557 weight_ents.reserve(rend - rstart + 1);
558
559 struct AdjBridge {
560 EntityHandle ent;
561 std::vector<int> adj;
562 AdjBridge(const EntityHandle ent, std::vector<int> &adj)
563 : ent(ent), adj(adj) {}
564 };
565
566 typedef multi_index_container<
567 AdjBridge,
568 indexed_by<
569
570 hashed_unique<member<AdjBridge, EntityHandle, &AdjBridge::ent>>
571
572 >>
573 AdjBridgeMap;
574
575 auto get_it = [&](auto i) {
576 auto it = ents.begin();
577 for (; i > 0; --i) {
578 if (it == ents.end())
579 break;
580 ++it;
581 }
582 return it;
583 };
584
585 Range proc_ents;
586 proc_ents.insert(get_it(rstart), get_it(rend));
587 if (proc_ents.size() != rend - rstart)
588 SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
589 "Wrong number of elements in range %d != %d", proc_ents.size(),
590 rend - rstart);
591
592 Range all_dim_ents;
593 CHKERR m_field.get_moab().get_adjacencies(
594 proc_ents, adj_dim, true, all_dim_ents, moab::Interface::UNION);
595
596 AdjBridgeMap adj_bridge_map;
597 auto hint = adj_bridge_map.begin();
598 std::vector<int> adj;
599 for (auto ent : all_dim_ents) {
600 Range adj_ents;
601 CHKERR m_field.get_moab().get_adjacencies(&ent, 1, dim, false, adj_ents);
602 adj_ents = intersect(adj_ents, ents);
603 adj.clear();
604 adj.reserve(adj_ents.size());
605 for (auto a : adj_ents)
606 adj.emplace_back(ents.index(a));
607 hint = adj_bridge_map.emplace_hint(hint, ent, adj);
608 }
609
610 int *_i;
611 int *_j;
612 {
613 const int nb_loc_elements = rend - rstart;
614 std::vector<int> i(nb_loc_elements + 1, 0), j;
615 {
616 std::vector<int> row_adj;
617 Range::iterator fe_it;
618 int ii, jj;
619 size_t max_row_size;
620 for (fe_it = proc_ents.begin(), ii = rstart, jj = 0, max_row_size = 0;
621 fe_it != proc_ents.end(); ++fe_it, ++ii) {
622
623 if (type_from_handle(*fe_it) == MBENTITYSET) {
624 SETERRQ(
625 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
626 "not yet implemented, don't know what to do for meshset element");
627 } else {
628
629 Range dim_ents;
630 CHKERR m_field.get_moab().get_adjacencies(&*fe_it, 1, adj_dim, false,
631 dim_ents);
632 dim_ents = intersect(dim_ents, all_dim_ents);
633
634 row_adj.clear();
635 for (auto e : dim_ents) {
636 auto adj_it = adj_bridge_map.find(e);
637 if (adj_it != adj_bridge_map.end()) {
638
639 for (const auto idx : adj_it->adj)
640 row_adj.push_back(idx);
641
642 } else
643 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
644 "Entity not found");
645 }
646
647 std::sort(row_adj.begin(), row_adj.end());
648 auto end = std::unique(row_adj.begin(), row_adj.end());
649
650 size_t row_size = std::distance(row_adj.begin(), end);
651 max_row_size = std::max(max_row_size, row_size);
652 if (j.capacity() < (nb_loc_elements - jj) * max_row_size)
653 j.reserve(nb_loc_elements * max_row_size);
654
655 i[jj] = j.size();
656 auto diag = ents.index(*fe_it);
657 for (auto it = row_adj.begin(); it != end; ++it)
658 if (*it != diag)
659 j.push_back(*it);
660 }
661
662 ++jj;
663
664 if (th_vertex_weights != NULL)
665 weight_ents.push_back(*fe_it);
666 }
667
668 i[jj] = j.size();
669 }
670
671 CHKERR PetscMalloc(i.size() * sizeof(int), &_i);
672 CHKERR PetscMalloc(j.size() * sizeof(int), &_j);
673 copy(i.begin(), i.end(), _i);
674 copy(j.begin(), j.end(), _j);
675 }
676
677 // get weights
678 int *vertex_weights = NULL;
679 if (th_vertex_weights != NULL) {
680 CHKERR PetscMalloc(weight_ents.size() * sizeof(int), &vertex_weights);
681 CHKERR m_field.get_moab().tag_get_data(*th_vertex_weights,
682 &*weight_ents.begin(),
683 weight_ents.size(), vertex_weights);
684 }
685
686 {
687 Mat Adj;
688 // Adjacency matrix used to partition problems, f.e. METIS
689 CHKERR MatCreateMPIAdj(m_field.get_comm(), rend - rstart, nb_elems, _i, _j,
690 PETSC_NULL, &Adj);
691 CHKERR MatSetOption(Adj, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
692
693 if (debug) {
694 Mat A;
695 CHKERR MatConvert(Adj, MATMPIAIJ, MAT_INITIAL_MATRIX, &A);
696 CHKERR MatView(A, PETSC_VIEWER_DRAW_WORLD);
697 std::string wait;
698 std::cin >> wait;
699 CHKERR MatDestroy(&A);
700 }
701
702 // run pets to do partitioning
703 MOFEM_LOG("WORLD", Sev::verbose) << "Start";
704
705 MatPartitioning part;
706 IS is;
707 CHKERR MatPartitioningCreate(m_field.get_comm(), &part);
708
709 CHKERR MatPartitioningSetAdjacency(part, Adj);
710 CHKERR MatPartitioningSetFromOptions(part);
711 CHKERR MatPartitioningSetNParts(part, n_parts);
712 if (th_vertex_weights != NULL) {
713 CHKERR MatPartitioningSetVertexWeights(part, vertex_weights);
714 }
715 PetscBool same;
716 PetscObjectTypeCompare((PetscObject)part, MATPARTITIONINGPARMETIS, &same);
717 if (same) {
718#ifdef PARMETIS
719 CHKERR MatPartitioningApply_Parmetis_MoFEM(part, &is);
720#endif
721 } else {
722 CHKERR MatPartitioningApply(part, &is);
723 }
724
725 MOFEM_LOG("WORLD", Sev::verbose) << "End";
726
727 // gather
728 IS is_gather, is_num, is_gather_num;
729 CHKERR ISAllGather(is, &is_gather);
730 CHKERR ISPartitioningToNumbering(is, &is_num);
731 CHKERR ISAllGather(is_num, &is_gather_num);
732
733 const int *part_number, *gids;
734 CHKERR ISGetIndices(is_gather, &part_number);
735 CHKERR ISGetIndices(is_gather_num, &gids);
736
737 // set partition tag and gid tag to entities
738 ParallelComm *pcomm = ParallelComm::get_pcomm(
739 &m_field.get_moab(), m_field.get_basic_entity_data_ptr()->pcommID);
740 Tag part_tag = pcomm->part_tag();
741 CHKERR m_field.get_moab().tag_set_data(part_tag, ents, part_number);
742 Tag gid_tag = m_field.get_moab().globalId_tag();
743
744 std::map<int, Range> parts_ents;
745 {
746 // get entities on each part
747 Range::iterator eit = ents.begin();
748 for (int ii = 0; eit != ents.end(); eit++, ii++) {
749 parts_ents[part_number[ii]].insert(*eit);
750 }
751 Range tagged_sets;
752 CHKERR m_field.get_moab().get_entities_by_type_and_tag(
753 0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
754 moab::Interface::UNION);
755 if (!tagged_sets.empty())
756 CHKERR m_field.get_moab().tag_delete_data(part_tag, tagged_sets);
757
758 if (n_parts > (int)tagged_sets.size()) {
759 // too few partition sets - create missing ones
760 int num_new = n_parts - tagged_sets.size();
761 for (int i = 0; i < num_new; i++) {
762 EntityHandle new_set;
763 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, new_set);
764 tagged_sets.insert(new_set);
765 }
766 } else if (n_parts < (int)tagged_sets.size()) {
767 // too many partition sets - delete extras
768 int num_del = tagged_sets.size() - n_parts;
769 for (int i = 0; i < num_del; i++) {
770 EntityHandle old_set = tagged_sets.pop_back();
771 CHKERR m_field.get_moab().delete_entities(&old_set, 1);
772 }
773 }
774 // write a tag to those sets denoting they're partition sets, with a value
775 // of the proc number
776 std::vector<int> dum_ids(n_parts);
777 for (int i = 0; i < n_parts; i++)
778 dum_ids[i] = i;
779 CHKERR m_field.get_moab().tag_set_data(part_tag, tagged_sets,
780 &*dum_ids.begin());
781 CHKERR m_field.get_moab().clear_meshset(tagged_sets);
782
783 // get lower dimension entities on each part
784 for (int pp = 0; pp != n_parts; pp++) {
785 Range dim_ents = parts_ents[pp].subset_by_dimension(dim);
786 for (int dd = dim - 1; dd >= 0; dd--) {
787 Range adj_ents;
788 if (dd > 0) {
789 CHKERR m_field.get_moab().get_adjacencies(
790 dim_ents, dd, false, adj_ents, moab::Interface::UNION);
791 } else {
792 CHKERR m_field.get_moab().get_connectivity(dim_ents, adj_ents,
793 true);
794 }
795 parts_ents[pp].merge(adj_ents);
796 }
797 }
798 for (int pp = 1; pp != n_parts; pp++) {
799 for (int ppp = 0; ppp != pp; ppp++) {
800 parts_ents[pp] = subtract(parts_ents[pp], parts_ents[ppp]);
801 }
802 }
803 if (debug) {
804 if (m_field.get_comm_rank() == 0) {
805 for (int rr = 0; rr != n_parts; rr++) {
806 ostringstream ss;
807 ss << "out_part_" << rr << ".vtk";
808 MOFEM_LOG("SELF", Sev::inform)
809 << "Save debug part mesh " << ss.str();
810 EntityHandle meshset;
811 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
812 CHKERR m_field.get_moab().add_entities(meshset, parts_ents[rr]);
813 CHKERR m_field.get_moab().write_file(ss.str().c_str(), "VTK", "",
814 &meshset, 1);
815 CHKERR m_field.get_moab().delete_entities(&meshset, 1);
816 }
817 }
818 }
819 for (int pp = 0; pp != n_parts; pp++) {
820 CHKERR m_field.get_moab().add_entities(tagged_sets[pp], parts_ents[pp]);
821 }
822
823 // set gid and part tag
824 for (EntityType t = MBVERTEX; t != MBENTITYSET; ++t) {
825
826 void *ptr;
827 int count;
828
829 int gid = 1; // moab indexing from 1a
830 for (int pp = 0; pp != n_parts; pp++) {
831 Range type_ents = parts_ents[pp].subset_by_type(t);
832 if (t != MBVERTEX) {
833 CHKERR m_field.get_moab().tag_clear_data(part_tag, type_ents, &pp);
834 }
835
836 auto eit = type_ents.begin();
837 for (; eit != type_ents.end();) {
838 CHKERR m_field.get_moab().tag_iterate(gid_tag, eit, type_ents.end(),
839 count, ptr);
840 auto gid_tag_ptr = static_cast<int *>(ptr);
841 for (; count > 0; --count) {
842 *gid_tag_ptr = gid;
843 ++eit;
844 ++gid;
845 ++gid_tag_ptr;
846 }
847 }
848 }
849 }
850 }
851
852 CHKERR ISRestoreIndices(is_gather, &part_number);
853 CHKERR ISRestoreIndices(is_gather_num, &gids);
854 CHKERR ISDestroy(&is_num);
855 CHKERR ISDestroy(&is_gather_num);
856 CHKERR ISDestroy(&is_gather);
857 CHKERR ISDestroy(&is);
858 CHKERR MatPartitioningDestroy(&part);
859 CHKERR MatDestroy(&Adj);
860 }
861
863}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
static const bool debug
const int dim
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
double A
constexpr double t
plate stiffness
Definition: plate.cpp:76
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 91 of file ProblemsManager.cpp.

94 {
95 return static_cast<MoFEM::Interface &>(cOre)
96 .getInterface<CommInterface>()
97 ->partitionMesh(ents, dim, adj_dim, n_parts, th_vertex_weights,
98 th_edge_weights, th_part_weights, verb, debug);
99}

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

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

◆ partitionSimpleProblem()

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

partition problem dofs

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

Definition at line 1541 of file ProblemsManager.cpp.

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

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

3029 {
3030 MoFEM::Interface &m_field = cOre;
3032
3033 auto bit_manager = m_field.getInterface<BitRefManager>();
3034
3035 Range ents;
3036 if (ents_ptr) {
3037 ents = *ents_ptr;
3038 CHKERR bit_manager->filterEntitiesByRefLevel(bit_ref_level, bit_ref_mask,
3039 ents, verb);
3040 } else {
3041 CHKERR bit_manager->getEntitiesByRefLevel(bit_ref_level, bit_ref_mask, ents,
3042 verb);
3043 }
3044
3045 CHKERR removeDofsOnEntities(problem_name, field_name, ents, lo_coeff,
3046 hi_coeff, lo_order, hi_order, verb, debug);
3047
3049}
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 2600 of file ProblemsManager.cpp.

2603 {
2604
2605 MoFEM::Interface &m_field = cOre;
2607
2608 const Problem *prb_ptr;
2609 CHKERR m_field.get_problem(problem_name, &prb_ptr);
2610
2611 decltype(prb_ptr->numeredRowDofsPtr) numered_dofs[2] = {
2612 prb_ptr->numeredRowDofsPtr, nullptr};
2613 if (prb_ptr->numeredRowDofsPtr != prb_ptr->numeredColDofsPtr)
2614 numered_dofs[1] = prb_ptr->numeredColDofsPtr;
2615
2616 int *nbdof_ptr[] = {&prb_ptr->nbDofsRow, &prb_ptr->nbDofsCol};
2617 int *local_nbdof_ptr[] = {&prb_ptr->nbLocDofsRow, &prb_ptr->nbLocDofsCol};
2618 int *ghost_nbdof_ptr[] = {&prb_ptr->nbGhostDofsRow, &prb_ptr->nbGhostDofsCol};
2619
2620 const int nb_init_row_dofs = prb_ptr->getNbDofsRow();
2621 const int nb_init_col_dofs = prb_ptr->getNbDofsCol();
2622 const int nb_init_loc_row_dofs = prb_ptr->getNbLocalDofsRow();
2623 const int nb_init_loc_col_dofs = prb_ptr->getNbLocalDofsCol();
2624 const int nb_init_ghost_row_dofs = prb_ptr->getNbGhostDofsRow();
2625 const int nb_init_ghost_col_dofs = prb_ptr->getNbGhostDofsCol();
2626
2627 for (int s = 0; s != 2; ++s)
2628 if (numered_dofs[s]) {
2629
2630 typedef multi_index_container<
2631
2632 NumeredDofEntity_multiIndex::iterator, indexed_by<sequenced<>>
2633
2634 >
2635 NumeredDofEntity_it_view_multiIndex;
2636
2637 const auto bit_number = m_field.get_field_bit_number(field_name);
2638 NumeredDofEntity_it_view_multiIndex dofs_it_view;
2639
2640 // Set -1 to global and local dofs indices
2641 for (auto pit = ents.const_pair_begin(); pit != ents.const_pair_end();
2642 ++pit) {
2643 auto lo = numered_dofs[s]->get<Unique_mi_tag>().lower_bound(
2644 DofEntity::getLoFieldEntityUId(bit_number, pit->first));
2645 auto hi = numered_dofs[s]->get<Unique_mi_tag>().upper_bound(
2646 DofEntity::getHiFieldEntityUId(bit_number, pit->second));
2647
2648 for (; lo != hi; ++lo)
2649 if ((*lo)->getDofCoeffIdx() >= lo_coeff &&
2650 (*lo)->getDofCoeffIdx() <= hi_coeff &&
2651 (*lo)->getDofOrder() >= lo_order &&
2652 (*lo)->getDofOrder() <= hi_order)
2653 dofs_it_view.emplace_back(numered_dofs[s]->project<0>(lo));
2654 }
2655
2656 if (verb > QUIET) {
2657 for (auto &dof : dofs_it_view)
2658 MOFEM_LOG("SYNC", Sev::noisy) << **dof;
2659
2661 }
2662
2663 // set negative index
2664 auto mod = NumeredDofEntity_part_and_all_indices_change(-1, -1, -1, -1);
2665 for (auto dit : dofs_it_view) {
2666 bool success = numered_dofs[s]->modify(dit, mod);
2667 if (!success)
2668 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2669 "can not set negative indices");
2670 }
2671
2672 // create weak view
2673 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_weak_view;
2674 dofs_weak_view.reserve(dofs_it_view.size());
2675 for (auto dit : dofs_it_view)
2676 dofs_weak_view.push_back(*dit);
2677
2678 if (verb >= NOISY)
2679 MOFEM_LOG_C("SYNC", Sev::noisy,
2680 "Number of DOFs in multi-index %d and to delete %d\n",
2681 numered_dofs[s]->size(), dofs_it_view.size());
2682
2683 // erase dofs from problem
2684 for (auto weak_dit : dofs_weak_view)
2685 if (auto dit = weak_dit.lock()) {
2686 numered_dofs[s]->erase(dit->getLocalUniqueId());
2687 }
2688
2689 if (verb >= NOISY)
2690 MOFEM_LOG_C("SYNC", Sev::noisy,
2691 "Number of DOFs in multi-index after delete %d\n",
2692 numered_dofs[s]->size());
2693
2694 // get current number of ghost dofs
2695 int nb_local_dofs = 0;
2696 int nb_ghost_dofs = 0;
2697 for (auto dit = numered_dofs[s]->get<PetscLocalIdx_mi_tag>().begin();
2698 dit != numered_dofs[s]->get<PetscLocalIdx_mi_tag>().end(); ++dit) {
2699 if ((*dit)->getPetscLocalDofIdx() >= 0 &&
2700 (*dit)->getPetscLocalDofIdx() < *(local_nbdof_ptr[s]))
2701 ++nb_local_dofs;
2702 else if ((*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s]))
2703 ++nb_ghost_dofs;
2704 else
2705 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2706 "Impossible case. You could run problem on no distributed "
2707 "mesh. That is not implemented. Dof local index is %d",
2708 (*dit)->getPetscLocalDofIdx());
2709 }
2710
2711 // get indices
2712 auto get_indices_by_tag = [&](auto tag, auto &indices, bool only_local) {
2713 const int nb_dofs = numered_dofs[s]->size();
2714 indices.clear();
2715 indices.reserve(nb_dofs);
2716 for (auto dit = numered_dofs[s]->get<decltype(tag)>().begin();
2717 dit != numered_dofs[s]->get<decltype(tag)>().end(); ++dit) {
2718 bool add = true;
2719 if (only_local) {
2720 if ((*dit)->getPetscLocalDofIdx() < 0 ||
2721 (*dit)->getPetscLocalDofIdx() >= *(local_nbdof_ptr[s])) {
2722 add = false;
2723 }
2724 }
2725 if (add)
2726 indices.push_back(decltype(tag)::get_index(dit));
2727 }
2728 };
2729
2730 auto get_indices_by_uid = [&](auto tag, auto &indices) {
2731 const int nb_dofs = numered_dofs[s]->size();
2732 indices.clear();
2733 indices.reserve(nb_dofs);
2734 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2735 ++dit)
2736 indices.push_back(decltype(tag)::get_index(dit));
2737 };
2738
2739 auto concatenate_dofs = [&](auto tag, auto &indices,
2740 const auto local_only) {
2742 get_indices_by_tag(tag, indices, local_only);
2743
2744 AO ao;
2745 if (local_only)
2746 CHKERR AOCreateMapping(m_field.get_comm(), indices.size(),
2747 &*indices.begin(), PETSC_NULL, &ao);
2748 else
2749 CHKERR AOCreateMapping(PETSC_COMM_SELF, indices.size(),
2750 &*indices.begin(), PETSC_NULL, &ao);
2751
2752 get_indices_by_uid(tag, indices);
2753 CHKERR AOApplicationToPetsc(ao, indices.size(), &*indices.begin());
2754 CHKERR AODestroy(&ao);
2756 };
2757
2758 // set indices index
2759 auto set_concatinated_indices = [&]() {
2760 std::vector<int> global_indices;
2761 std::vector<int> local_indices;
2763 CHKERR concatenate_dofs(PetscGlobalIdx_mi_tag(), global_indices, true);
2764 CHKERR concatenate_dofs(PetscLocalIdx_mi_tag(), local_indices, false);
2765 auto gi = global_indices.begin();
2766 auto li = local_indices.begin();
2767 for (auto dit = numered_dofs[s]->begin(); dit != numered_dofs[s]->end();
2768 ++dit) {
2769 auto mod = NumeredDofEntity_part_and_all_indices_change(
2770 (*dit)->getPart(), (*dit)->getDofIdx(), *gi, *li);
2771 bool success = numered_dofs[s]->modify(dit, mod);
2772 if (!success)
2773 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2774 "can not set negative indices");
2775 ++gi;
2776 ++li;
2777 }
2779 };
2780 CHKERR set_concatinated_indices();
2781
2782 MPI_Allreduce(&nb_local_dofs, nbdof_ptr[s], 1, MPI_INT, MPI_SUM,
2783 m_field.get_comm());
2784 *(local_nbdof_ptr[s]) = nb_local_dofs;
2785 *(ghost_nbdof_ptr[s]) = nb_ghost_dofs;
2786
2787 if (debug)
2788 for (auto dof : (*numered_dofs[s])) {
2789 if (dof->getPetscGlobalDofIdx() < 0) {
2790 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2791 "Negative global idx");
2792 }
2793 if (dof->getPetscLocalDofIdx() < 0) {
2794 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
2795 "Negative local idx");
2796 }
2797 }
2798
2799 } else {
2800
2801 *(nbdof_ptr[1]) = *(nbdof_ptr[0]);
2802 *(local_nbdof_ptr[1]) = *(local_nbdof_ptr[0]);
2803 *(ghost_nbdof_ptr[1]) = *(ghost_nbdof_ptr[0]);
2804 }
2805
2806 if (verb > QUIET) {
2807 MOFEM_LOG_C("SYNC", Sev::verbose,
2808 "removed ents on rank %d from problem %s dofs [ %d / %d "
2809 "(before %d / "
2810 "%d) local, %d / %d (before %d / %d) "
2811 "ghost, %d / %d (before %d / %d) global]",
2812 m_field.get_comm_rank(), prb_ptr->getName().c_str(),
2813 prb_ptr->getNbLocalDofsRow(), prb_ptr->getNbLocalDofsCol(),
2814 nb_init_row_dofs, nb_init_col_dofs,
2815 prb_ptr->getNbGhostDofsRow(), prb_ptr->getNbGhostDofsCol(),
2816 nb_init_ghost_row_dofs, nb_init_ghost_col_dofs,
2817 prb_ptr->getNbDofsRow(), prb_ptr->getNbDofsCol(),
2818 nb_init_loc_row_dofs, nb_init_loc_col_dofs);
2820 }
2821
2823}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem