v0.13.1
Loading...
Searching...
No Matches
Files | Classes | Functions
ProblemsManager

Adding and managing problems. More...

Collaboration diagram for ProblemsManager:

Files

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

Classes

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

Functions

DEPRECATED MoFEMErrorCode MoFEM::ProblemsManager::partitionMesh (const Range &ents, const int dim, const int adj_dim, const int n_parts, Tag *th_vertex_weights=nullptr, Tag *th_edge_weights=nullptr, Tag *th_part_weights=nullptr, int verb=VERBOSE, const bool debug=false)
 Set partition tag to each finite element in the problem. More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblem (const std::string name, const bool square_matrix, int verb=VERBOSE)
 build problem data structures More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblem (Problem *problem_ptr, const bool square_matrix, int verb=VERBOSE)
 build problem data structures More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh (const std::string name, const bool square_matrix, int verb=VERBOSE)
 build problem data structures, assuming that mesh is distributed (collective) More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildProblemOnDistributedMesh (Problem *problem_ptr, const bool square_matrix=true, int verb=VERBOSE)
 build problem data structures, assuming that mesh is distributed (collective) More...
 
MoFEMErrorCode MoFEM::ProblemsManager::buildSubProblem (const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, 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 1238 of file ProblemsManager.cpp.

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

◆ buildProblem() [1/2]

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

build problem data structures

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

Definition at line 87 of file ProblemsManager.cpp.

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

◆ buildProblem() [2/2]

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

build problem data structures

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

Definition at line 107 of file ProblemsManager.cpp.

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

◆ buildProblemOnDistributedMesh() [1/2]

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

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

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

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

Examples
nonlinear_dynamics.cpp, and remove_entities_from_problem.cpp.

Definition at line 290 of file ProblemsManager.cpp.

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

◆ buildProblemOnDistributedMesh() [2/2]

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

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

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

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

Definition at line 313 of file ProblemsManager.cpp.

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

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

2551 {
2552 MoFEM::Interface &m_field = cOre;
2553 const Problem *problem_ptr;
2555
2556 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, *meshset);
2557 CHKERR m_field.get_problem(prb_name, &problem_ptr);
2558 auto fit =
2559 problem_ptr->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
2560 .lower_bound(fe_name);
2561 auto hi_fe_it =
2562 problem_ptr->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
2563 .upper_bound(fe_name);
2564 std::vector<EntityHandle> fe_vec;
2565 fe_vec.reserve(std::distance(fit, hi_fe_it));
2566 for (; fit != hi_fe_it; fit++)
2567 fe_vec.push_back(fit->get()->getEnt());
2568 CHKERR m_field.get_moab().add_entities(*meshset, &*fe_vec.begin(),
2569 fe_vec.size());
2571}
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 2574 of file ProblemsManager.cpp.

2576 {
2577 MoFEM::Interface &m_field = cOre;
2578 const Problem *problem_ptr;
2580 CHKERR m_field.get_problem(name, &problem_ptr);
2581 CHKERR problem_ptr->getNumberOfElementsByNameAndPart(PETSC_COMM_WORLD,
2582 fe_name, layout);
2584}

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

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

◆ partitionFiniteElements()

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

partition finite elements

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

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

Definition at line 2157 of file ProblemsManager.cpp.

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

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

◆ partitionGhostDofsOnDistributedMesh()

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

determine ghost nodes on distributed meshes

Parameters
nameproblem name

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

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

Definition at line 2463 of file ProblemsManager.cpp.

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

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

◆ partitionMesh() [2/2]

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

Set partition tag to each finite element in the problem.

Deprecated:
Use CommInterface

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

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

Definition at line 77 of file ProblemsManager.cpp.

80 {
81 return static_cast<MoFEM::Interface &>(cOre)
82 .getInterface<CommInterface>()
83 ->partitionMesh(ents, dim, adj_dim, n_parts, th_vertex_weights,
84 th_edge_weights, th_part_weights, verb, debug);
85}

◆ partitionProblem()

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

partition problem dofs (collective)

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

Definition at line 1673 of file ProblemsManager.cpp.

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

◆ partitionSimpleProblem()

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

partition problem dofs

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

Definition at line 1527 of file ProblemsManager.cpp.

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

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

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

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