71 unsigned nsub_mesh = sub_mesh_pt.size();
74 unsigned long n_element = 0;
75 unsigned long n_node = 0;
80 for (
unsigned imesh = 0; imesh < nsub_mesh; imesh++)
82 n_element += sub_mesh_pt[imesh]->nelement();
83 n_node += sub_mesh_pt[imesh]->nnode();
84 n_bound += sub_mesh_pt[imesh]->nboundary();
100 std::set<GeneralisedElement*> element_set_pt;
101 std::set<Node*> node_set_pt;
104 unsigned ibound_global = 0;
106 for (
unsigned imesh = 0; imesh < nsub_mesh; imesh++)
110 unsigned nel_before = 0;
111 unsigned long n_element = sub_mesh_pt[imesh]->nelement();
112 for (
unsigned long e = 0;
e < n_element;
e++)
115 element_set_pt.insert(el_pt);
117 unsigned nel_now = element_set_pt.size();
118 if (nel_now == nel_before)
120 std::ostringstream warning_stream;
121 warning_stream <<
"WARNING: " << std::endl
122 <<
"Element " <<
e <<
" in submesh " << imesh
123 <<
" is a duplicate \n and was ignored when assembling"
124 <<
" combined mesh." << std::endl;
126 "Mesh::Mesh(const Vector<Mesh*>&)",
127 OOMPH_EXCEPTION_LOCATION);
133 nel_before = nel_now;
138 unsigned nnod_before = 0;
139 unsigned long n_node = sub_mesh_pt[imesh]->nnode();
140 for (
unsigned long n = 0; n < n_node; n++)
142 Node* nod_pt = sub_mesh_pt[imesh]->node_pt(n);
143 node_set_pt.insert(nod_pt);
145 unsigned nnod_now = node_set_pt.size();
146 if (nnod_now == nnod_before)
148 std::ostringstream warning_stream;
150 <<
"WARNING: " << std::endl
151 <<
"Node " << n <<
" in submesh " << imesh
152 <<
" is a duplicate \n and was ignored when assembling "
153 <<
"combined mesh." << std::endl;
155 "Mesh::Mesh(const Vector<Mesh*>&)",
156 OOMPH_EXCEPTION_LOCATION);
162 nnod_before = nnod_now;
166 unsigned n_bound = sub_mesh_pt[imesh]->nboundary();
167 for (
unsigned ibound = 0; ibound < n_bound; ibound++)
171 unsigned long n_bound_node = sub_mesh_pt[imesh]->nboundary_node(ibound);
172 for (
unsigned long n = 0; n < n_bound_node; n++)
193 for (
unsigned n = 0; n < n_boundary_node; n++)
208 for (
unsigned b = 0; b < n_bound; b++)
252 bool node_already_on_this_boundary =
false;
254 for (
unsigned n = 0; n < nbound_node; n++)
259 node_already_on_this_boundary =
true;
264 if (!node_already_on_this_boundary)
296 for (
unsigned long n = 0; n <
nnode(); n++)
303 for (
unsigned imaster = 0; imaster < nmaster; imaster++)
325 std::ostringstream err_stream;
327 err_stream <<
"Calling node_update() for a mesh which contains"
329 <<
"master nodes which live in the external storage."
331 <<
"These nodes do not belong to elements on this"
333 <<
"processor and therefore cannot be updated locally."
337 OOMPH_CURRENT_FUNCTION,
338 OOMPH_EXCEPTION_LOCATION);
356 std::map<Node*, bool> has_node_been_updated;
359 unsigned n_node =
nnode();
362 for (
unsigned n = 0; n < n_node; n++)
368 has_node_been_updated[nod_pt] =
false;
373 for (
unsigned e = 0;
e < nel;
e++)
383 unsigned ndim_el = el_pt->
dim();
387 unsigned n_node = el_pt->
nnode();
388 for (
unsigned j = 0; j < n_node; j++)
394 unsigned ndim_node = nod_pt->
ndim();
401 if (!has_node_been_updated[nod_pt])
413 for (
unsigned i = 0;
i < ndim_node;
i++)
416 if (solid_node_pt != 0)
419 if (update_all_solid_nodes)
421 solid_node_pt->
x(
i) = r[
i];
433 has_node_been_updated[nod_pt] =
true;
452 unsigned nnod = ext_halo_node_pt.size();
453 for (
unsigned j = 0; j < nnod; j++)
455 ext_halo_node_pt[j]->node_update();
462 for (
unsigned long n = 0; n < n_node; n++)
468 unsigned ndim_node = nod_pt->
ndim();
471 for (
unsigned i = 0;
i < ndim_node;
i++)
478 for (
unsigned imaster = 0; imaster < nmaster; imaster++)
481 for (
unsigned i = 0;
i < ndim_node;
i++)
493 for (
unsigned long n = 0; n < n_node; n++)
495 Node_pt[n]->perform_auxiliary_node_update_fct();
499 oomph_info <<
"Time taken to update nodal positions [sec]: "
517 unsigned n_node =
nnode();
520 for (
unsigned i = 0;
i < n_node;
i++)
533 const bool& use_old_ordering)
const
536 if (use_old_ordering)
539 std::map<Node*, bool> done;
542 unsigned nnod =
nnode();
545 reordering.assign(nnod, 0);
558 for (
unsigned j = 0; j < nnod; j++)
565 unsigned long count = 0;
571 for (
unsigned e = 0;
e < nel;
e++)
576 checked_dynamic_cast<FiniteElement*>(
element_pt(
e));
579 unsigned nnod = el_pt->
nnode();
582 for (
unsigned j = 0; j < nnod; j++)
599 reordering[count] = nod_pt;
614 std::string error_message =
"Trouble: Number of nodes hasn't stayed ";
617 error_message +=
"constant during reordering!\n";
621 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
627 unsigned n_node =
nnode();
630 reordering.resize(n_node);
633 for (
unsigned i = 0;
i < n_node;
i++)
640 std::sort(reordering.begin(),
654 unsigned long Node_pt_range =
Node_pt.size();
655 for (
unsigned long i = Node_pt_range;
i > 0;
i--)
663 unsigned long Element_pt_range =
Element_pt.size();
664 for (
unsigned long i = Element_pt_range;
i > 0;
i--)
680 unsigned long equation_number = Dof_pt.size();
683 unsigned long nnod =
Node_pt.size();
685 for (
unsigned long i = 0;
i < nnod;
i++)
687 Node_pt[
i]->assign_eqn_numbers(equation_number, Dof_pt);
692 for (
unsigned long i = 0;
i < nel;
i++)
694 Element_pt[
i]->assign_internal_eqn_numbers(equation_number, Dof_pt);
698 return (equation_number);
715 unsigned long nnod =
Node_pt.size();
716 for (
unsigned long i = 0;
i < nnod;
i++)
718 std::stringstream conversion;
719 conversion <<
" of Node " <<
i << current_string;
726 for (
unsigned long i = 0;
i < nel;
i++)
728 std::stringstream conversion;
729 conversion <<
" in Element " <<
i <<
" [" <<
typeid(*
Element_pt[
i]).name()
730 <<
"] " << current_string;
751 for (
unsigned long i = 0;
i < nel;
i++)
753 std::stringstream conversion;
754 conversion <<
" in Element" <<
i <<
" [" <<
typeid(*
Element_pt[
i]).name()
755 <<
"] " << current_string;
768 unsigned long Element_pt_range =
Element_pt.size();
769 for (
unsigned long i = 0;
i < Element_pt_range;
i++)
771 Element_pt[
i]->assign_local_eqn_numbers(store_local_dof_pt);
801 std::set<GeneralisedElement*> element_set_pt;
802 unsigned long Element_pt_range =
Element_pt.size();
803 for (
unsigned long i = 0;
i < Element_pt_range;
i++)
808 oomph_info <<
"\n ERROR: Failed Element::self_test() for element i="
816 if (element_set_pt.size() != Element_pt_range)
818 oomph_info <<
"ERROR: " << Element_pt_range - element_set_pt.size()
819 <<
" duplicate elements were encountered in mesh!"
826 std::set<Node*> node_set_pt;
827 unsigned long Node_pt_range =
Node_pt.size();
828 for (
unsigned long i = 0;
i < Node_pt_range;
i++)
833 oomph_info <<
"\n ERROR: Failed Node::self_test() for node i=" <<
i
841 if (node_set_pt.size() != Node_pt_range)
843 oomph_info <<
"ERROR: " << Node_pt_range - node_set_pt.size()
844 <<
" duplicate nodes were encountered in mesh!" << std::endl;
871 std::ofstream& inverted_element_file)
874 mesh_has_inverted_elements =
false;
883 for (
unsigned e = 0;
e < nelem;
e++)
891 unsigned n_node = el_pt->
nnode();
892 unsigned n_dim = el_pt->
dim();
897 if (n_dim == ndim_node)
901 DShape dpsidx(n_node, n_dim);
905 bool is_inverted =
false;
908 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
911 for (
unsigned i = 0;
i < n_dim;
i++)
946 mesh_has_inverted_elements =
true;
947 if (inverted_element_file.is_open())
949 el_pt->
output(inverted_element_file);
976 unsigned long n_node =
nnode();
977 for (
unsigned long n = 0; n < n_node; n++)
980 if (!(
Node_pt[n]->is_obsolete()))
982 new_node_pt.push_back(
Node_pt[n]);
990 if (!(
Node_pt[n]->is_on_boundary()))
992 deleted_node_pt.push_back(
Node_pt[n]);
1010 for (
unsigned ibound = 0; ibound < num_bound; ibound++)
1020 new_boundary_node_pt.reserve(Nboundary_node);
1022 for (
unsigned long n = 0; n < Nboundary_node; n++)
1039 deleted_node_pt.push_back(
1053 return deleted_node_pt;
1068 for (
unsigned long ibound = 0; ibound < num_bound; ibound++)
1073 outfile <<
"ZONE T=\"boundary" << ibound <<
"\"\n";
1075 for (
unsigned inod = 0; inod < nnod; inod++)
1088 void Mesh::dump(std::ofstream& dump_file,
const bool& use_old_ordering)
const
1096 unsigned long Node_pt_range = this->
nnode();
1099 dump_file << Node_pt_range <<
" # number of nodes " << std::endl;
1102 for (
unsigned nd = 0; nd < Node_pt_range; nd++)
1104 reordering[nd]->dump(dump_file);
1108 unsigned n_element = this->
nelement();
1109 for (
unsigned e = 0;
e < n_element;
e++)
1115 dump_file << n_internal
1116 <<
" # number of internal Data items in element " <<
e
1118 for (
unsigned i = 0;
i < n_internal;
i++)
1142 unsigned long n_node = this->
nnode();
1145 getline(restart_file, input_string,
'#');
1148 restart_file.ignore(80,
'\n');
1151 unsigned long check_n_node = atoi(input_string.c_str());
1152 if (check_n_node != n_node)
1154 std::ostringstream error_stream;
1155 error_stream <<
"The number of nodes allocated " << n_node
1156 <<
" is not the same as specified in the restart file "
1157 << check_n_node << std::endl;
1160 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1164 for (
unsigned long n = 0; n < n_node; n++)
1168 if (el_node_pt != 0)
1170 el_node_pt->
read(restart_file);
1181 unsigned n_element = this->
nelement();
1182 for (
unsigned e = 0;
e < n_element;
e++)
1189 getline(restart_file, input_string,
'#');
1192 restart_file.ignore(80,
'\n');
1195 unsigned long check_n_internal = atoi(input_string.c_str());
1196 if (check_n_internal != n_internal)
1198 std::ostringstream error_stream;
1199 error_stream <<
"The number of internal data " << n_internal
1200 <<
" is not the same as specified in the restart file "
1201 << check_n_internal << std::endl;
1204 OOMPH_CURRENT_FUNCTION,
1205 OOMPH_EXCEPTION_LOCATION);
1208 for (
unsigned i = 0;
i < n_internal;
i++)
1226 const unsigned& nplot)
const
1229 file_out.setf(std::ios_base::uppercase);
1232 unsigned long number_of_elements = this->
Element_pt.size();
1240 throw OomphLibError(
"Recast for FiniteElement failed for element 0!\n",
1241 OOMPH_CURRENT_FUNCTION,
1242 OOMPH_EXCEPTION_LOCATION);
1251 for (
unsigned i = 1;
i < number_of_elements;
i++)
1255 if (el_zero_ndof != el_i_ndof)
1257 std::stringstream error_stream;
1259 <<
"Element " <<
i <<
" has different number of degrees of freedom\n"
1260 <<
"than from previous elements, Paraview cannot handle this.\n"
1261 <<
"We suggest that the problem is broken up into submeshes instead."
1264 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1270 unsigned long number_of_nodes = 0;
1271 unsigned long total_number_of_elements = 0;
1274 for (
unsigned i = 0;
i < number_of_elements;
i++)
1283 std::stringstream error_stream;
1284 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
1286 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1300 file_out <<
"<?xml version=\"1.0\"?>\n"
1301 <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1302 <<
"byte_order=\"LittleEndian\">\n"
1303 <<
"<UnstructuredGrid>\n"
1304 <<
"<Piece NumberOfPoints=\"" << number_of_nodes
1305 <<
"\" NumberOfCells=\"" << total_number_of_elements <<
"\">\n";
1315 file_out <<
"<PointData ";
1323 for (
unsigned i = 0;
i < ndof;
i++)
1325 file_out <<
"<DataArray type=\"Float32\" "
1327 <<
"format=\"ascii\""
1330 for (
unsigned j = 0; j < number_of_elements; j++)
1339 std::stringstream error_stream;
1340 error_stream <<
"Recast for element " << j <<
" failed" << std::endl;
1342 OOMPH_CURRENT_FUNCTION,
1343 OOMPH_EXCEPTION_LOCATION);
1351 file_out <<
"</DataArray>\n";
1355 file_out <<
"</PointData>\n";
1361 file_out <<
"<Points>\n"
1362 <<
"<DataArray type=\"Float32\""
1363 <<
" NumberOfComponents=\""
1366 <<
"format=\"ascii\">\n";
1369 for (
unsigned i = 0;
i < number_of_elements;
i++)
1378 std::stringstream error_stream;
1379 error_stream <<
"Recast for element " <<
i <<
" faild" << std::endl;
1381 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1388 file_out <<
"</DataArray>\n"
1397 <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
1401 unsigned counter = 0;
1404 for (
unsigned i = 0;
i < number_of_elements;
i++)
1413 std::stringstream error_stream;
1414 error_stream <<
"Recast for element " <<
i <<
" faild" << std::endl;
1416 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1422 file_out <<
"</DataArray>\n"
1423 <<
"<DataArray type=\"Int32\" "
1424 <<
"Name=\"offsets\" format=\"ascii\">\n";
1427 unsigned offset_sum = 0;
1430 for (
unsigned i = 0;
i < number_of_elements;
i++)
1439 std::stringstream error_stream;
1440 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
1442 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1448 file_out <<
"</DataArray>\n"
1449 <<
"<DataArray type=\"UInt8\" Name=\"types\">\n";
1452 for (
unsigned i = 0;
i < number_of_elements;
i++)
1461 std::stringstream error_stream;
1462 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
1464 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1471 file_out <<
"</DataArray>\n"
1477 file_out <<
"</Piece>\n"
1478 <<
"</UnstructuredGrid>\n"
1492 std::ofstream& file_out,
1493 const unsigned& nplot,
1497 file_out.setf(std::ios_base::uppercase);
1500 unsigned long number_of_elements = this->
Element_pt.size();
1508 throw OomphLibError(
"Recast for FiniteElement failed for element 0!\n",
1509 OOMPH_CURRENT_FUNCTION,
1510 OOMPH_EXCEPTION_LOCATION);
1519 for (
unsigned i = 1;
i < number_of_elements;
i++)
1523 if (el_zero_ndof != el_i_ndof)
1525 std::stringstream error_stream;
1527 <<
"Element " <<
i <<
" has different number of degrees of freedom\n"
1528 <<
"than from previous elements, Paraview cannot handle this.\n"
1529 <<
"We suggest that the problem is broken up into submeshes instead."
1532 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1538 unsigned long number_of_nodes = 0;
1539 unsigned long total_number_of_elements = 0;
1542 for (
unsigned i = 0;
i < number_of_elements;
i++)
1551 std::stringstream error_stream;
1552 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
1554 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1568 file_out <<
"<?xml version=\"1.0\"?>\n"
1569 <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1570 <<
"byte_order=\"LittleEndian\">\n"
1571 <<
"<UnstructuredGrid>\n"
1572 <<
"<Piece NumberOfPoints=\"" << number_of_nodes
1573 <<
"\" NumberOfCells=\"" << total_number_of_elements <<
"\">\n";
1583 file_out <<
"<PointData ";
1591 for (
unsigned i = 0;
i < ndof;
i++)
1593 file_out <<
"<DataArray type=\"Float32\" "
1595 <<
"format=\"ascii\""
1598 for (
unsigned j = 0; j < number_of_elements; j++)
1607 std::stringstream error_stream;
1608 error_stream <<
"Recast for element " << j <<
" failed" << std::endl;
1610 OOMPH_CURRENT_FUNCTION,
1611 OOMPH_EXCEPTION_LOCATION);
1619 file_out <<
"</DataArray>\n";
1623 file_out <<
"</PointData>\n";
1629 file_out <<
"<Points>\n"
1630 <<
"<DataArray type=\"Float32\""
1631 <<
" NumberOfComponents=\""
1634 <<
"format=\"ascii\">\n";
1637 for (
unsigned i = 0;
i < number_of_elements;
i++)
1646 std::stringstream error_stream;
1647 error_stream <<
"Recast for element " <<
i <<
" faild" << std::endl;
1649 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1656 file_out <<
"</DataArray>\n"
1665 <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
1669 unsigned counter = 0;
1672 for (
unsigned i = 0;
i < number_of_elements;
i++)
1681 std::stringstream error_stream;
1682 error_stream <<
"Recast for element " <<
i <<
" faild" << std::endl;
1684 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1690 file_out <<
"</DataArray>\n"
1691 <<
"<DataArray type=\"Int32\" "
1692 <<
"Name=\"offsets\" format=\"ascii\">\n";
1695 unsigned offset_sum = 0;
1698 for (
unsigned i = 0;
i < number_of_elements;
i++)
1707 std::stringstream error_stream;
1708 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
1710 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1716 file_out <<
"</DataArray>\n"
1717 <<
"<DataArray type=\"UInt8\" Name=\"types\">\n";
1720 for (
unsigned i = 0;
i < number_of_elements;
i++)
1729 std::stringstream error_stream;
1730 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
1732 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1739 file_out <<
"</DataArray>\n"
1745 file_out <<
"</Piece>\n"
1746 <<
"</UnstructuredGrid>\n"
1760 std::ofstream& file_out,
1761 const unsigned& nplot,
1766 file_out.setf(std::ios_base::uppercase);
1769 unsigned long number_of_elements = this->
Element_pt.size();
1777 throw OomphLibError(
"Recast for FiniteElement failed for element 0!\n",
1778 OOMPH_CURRENT_FUNCTION,
1779 OOMPH_EXCEPTION_LOCATION);
1788 for (
unsigned i = 1;
i < number_of_elements;
i++)
1792 if (el_zero_ndof != el_i_ndof)
1794 std::stringstream error_stream;
1796 <<
"Element " <<
i <<
" has different number of degrees of freedom\n"
1797 <<
"than from previous elements, Paraview cannot handle this.\n"
1798 <<
"We suggest that the problem is broken up into submeshes instead."
1801 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1807 unsigned long number_of_nodes = 0;
1808 unsigned long total_number_of_elements = 0;
1811 for (
unsigned i = 0;
i < number_of_elements;
i++)
1820 std::stringstream error_stream;
1821 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
1823 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1837 file_out <<
"<?xml version=\"1.0\"?>\n"
1838 <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1839 <<
"byte_order=\"LittleEndian\">\n"
1840 <<
"<UnstructuredGrid>\n"
1841 <<
"<Piece NumberOfPoints=\"" << number_of_nodes
1842 <<
"\" NumberOfCells=\"" << total_number_of_elements <<
"\">\n";
1852 file_out <<
"<PointData ";
1860 for (
unsigned i = 0;
i < ndof;
i++)
1862 file_out <<
"<DataArray type=\"Float32\" "
1864 <<
"format=\"ascii\""
1867 for (
unsigned j = 0; j < number_of_elements; j++)
1876 std::stringstream error_stream;
1877 error_stream <<
"Recast for element " << j <<
" failed" << std::endl;
1879 OOMPH_CURRENT_FUNCTION,
1880 OOMPH_EXCEPTION_LOCATION);
1885 file_out,
i, nplot, time, exact_soln_pt);
1889 file_out <<
"</DataArray>\n";
1893 file_out <<
"</PointData>\n";
1899 file_out <<
"<Points>\n"
1900 <<
"<DataArray type=\"Float32\""
1901 <<
" NumberOfComponents=\""
1904 <<
"format=\"ascii\">\n";
1907 for (
unsigned i = 0;
i < number_of_elements;
i++)
1916 std::stringstream error_stream;
1917 error_stream <<
"Recast for element " <<
i <<
" faild" << std::endl;
1919 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1926 file_out <<
"</DataArray>\n"
1935 <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
1939 unsigned counter = 0;
1942 for (
unsigned i = 0;
i < number_of_elements;
i++)
1951 std::stringstream error_stream;
1952 error_stream <<
"Recast for element " <<
i <<
" faild" << std::endl;
1954 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1960 file_out <<
"</DataArray>\n"
1961 <<
"<DataArray type=\"Int32\" "
1962 <<
"Name=\"offsets\" format=\"ascii\">\n";
1965 unsigned offset_sum = 0;
1968 for (
unsigned i = 0;
i < number_of_elements;
i++)
1977 std::stringstream error_stream;
1978 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
1980 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1986 file_out <<
"</DataArray>\n"
1987 <<
"<DataArray type=\"UInt8\" Name=\"types\">\n";
1990 for (
unsigned i = 0;
i < number_of_elements;
i++)
1999 std::stringstream error_stream;
2000 error_stream <<
"Recast for element " <<
i <<
" failed" << std::endl;
2002 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2009 file_out <<
"</DataArray>\n"
2015 file_out <<
"</Piece>\n"
2016 <<
"</UnstructuredGrid>\n"
2031 unsigned long Element_pt_range =
Element_pt.size();
2032 for (
unsigned long e = 0;
e < Element_pt_range;
e++)
2038 oomph_info <<
"Can't execute output(...) for non FiniteElements"
2043 #ifdef OOMPH_HAS_MPI
2049 #ifdef OOMPH_HAS_MPI
2073 unsigned long Element_pt_range =
Element_pt.size();
2075 for (
unsigned long e = 0;
e < Element_pt_range;
e++)
2081 oomph_info <<
"Can't execute output(...) for non FiniteElements"
2086 #ifdef OOMPH_HAS_MPI
2090 el_pt->
output(outfile, n_plot);
2092 #ifdef OOMPH_HAS_MPI
2097 el_pt->
output(outfile, n_plot);
2117 unsigned long Element_pt_range =
Element_pt.size();
2118 for (
unsigned long e = 0;
e < Element_pt_range;
e++)
2124 oomph_info <<
"Can't execute output(...) for non FiniteElements"
2129 #ifdef OOMPH_HAS_MPI
2135 #ifdef OOMPH_HAS_MPI
2160 unsigned long Element_pt_range =
Element_pt.size();
2161 for (
unsigned long e = 0;
e < Element_pt_range;
e++)
2167 oomph_info <<
"Can't execute output(...) for non FiniteElements"
2172 #ifdef OOMPH_HAS_MPI
2176 el_pt->
output(file_pt, n_plot);
2178 #ifdef OOMPH_HAS_MPI
2183 el_pt->
output(file_pt, n_plot);
2200 const unsigned& n_plot,
2205 unsigned long Element_pt_range =
Element_pt.size();
2206 for (
unsigned long e = 0;
e < Element_pt_range;
e++)
2212 oomph_info <<
"Can't execute output_fct(...) for non FiniteElements"
2217 #ifdef OOMPH_HAS_MPI
2221 el_pt->
output_fct(outfile, n_plot, exact_soln_pt);
2223 #ifdef OOMPH_HAS_MPI
2228 el_pt->
output_fct(outfile, n_plot, exact_soln_pt);
2244 const unsigned& n_plot,
2250 unsigned long Element_pt_range =
Element_pt.size();
2251 for (
unsigned long e = 0;
e < Element_pt_range;
e++)
2257 oomph_info <<
"Can't execute output_fct(...) for non FiniteElements"
2262 #ifdef OOMPH_HAS_MPI
2266 el_pt->
output_fct(outfile, n_plot, time, exact_soln_pt);
2268 #ifdef OOMPH_HAS_MPI
2273 el_pt->
output_fct(outfile, n_plot, time, exact_soln_pt);
2291 unsigned long Nelement =
nelement();
2292 for (
unsigned long e = 0;
e < Nelement;
e++)
2295 unsigned Ninternal =
element_pt(
e)->ninternal_data();
2298 for (
unsigned j = 0; j < Ninternal; j++)
2301 ->internal_data_pt(j)
2303 ->assign_initial_values_impulsive(
element_pt(
e)->internal_data_pt(j));
2308 unsigned long n_node =
nnode();
2309 for (
unsigned long n = 0; n < n_node; n++)
2312 Node_pt[n]->time_stepper_pt()->assign_initial_values_impulsive(
2316 ->position_time_stepper_pt()
2317 ->assign_initial_positions_impulsive(
Node_pt[n]);
2330 const unsigned long Nelement =
nelement();
2331 for (
unsigned long e = 0;
e < Nelement;
e++)
2334 const unsigned Ninternal =
element_pt(
e)->ninternal_data();
2337 for (
unsigned j = 0; j < Ninternal; j++)
2340 ->internal_data_pt(j)
2342 ->shift_time_values(
element_pt(
e)->internal_data_pt(j));
2347 const unsigned long n_node =
nnode();
2348 for (
unsigned long n = 0; n < n_node; n++)
2354 Node_pt[n]->position_time_stepper_pt()->shift_time_positions(
Node_pt[n]);
2370 const unsigned long Nelement =
nelement();
2371 for (
unsigned long e = 0;
e < Nelement;
e++)
2374 const unsigned Ninternal =
element_pt(
e)->ninternal_data();
2377 for (
unsigned j = 0; j < Ninternal; j++)
2380 ->internal_data_pt(j)
2382 ->calculate_predicted_values(
element_pt(
e)->internal_data_pt(j));
2387 const unsigned long n_node =
nnode();
2388 for (
unsigned long n = 0; n < n_node; n++)
2391 Node_pt[n]->time_stepper_pt()->calculate_predicted_values(
Node_pt[n]);
2393 Node_pt[n]->position_time_stepper_pt()->calculate_predicted_positions(
2403 const bool& preserve_existing_data)
2408 std::ostringstream warning_stream;
2410 <<
"Empty set_mesh_level_time_stepper() has been called.\n"
2411 <<
"This function needs to be overloaded to reset any (pointers to) \n"
2412 <<
"timesteppers for meshes that store timesteppers in locations "
2414 <<
"than the Nodes or Elements;\n"
2415 <<
"e.g. SpineMeshes have SpineData with timesteppers,\n"
2416 <<
"Triangle and TetMeshes store the timestepper for use in "
2417 "adaptivity.\n\n\n";
2419 <<
"If you are solving a continuation or bifurcation detecion\n"
2420 <<
"problem and strange things are happening, then check that\n"
2421 <<
"you don't need to overload this function for your mesh."
2422 <<
"\n This warning can be suppressed by setting:\n"
2423 <<
"Mesh::Suppress_warning_about_empty_mesh_level_time_stepper_"
2427 warning_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2440 const unsigned long n_node = this->
nnode();
2441 for (
unsigned long n = 0; n < n_node; n++)
2449 const unsigned long n_element = this->
nelement();
2450 for (
unsigned long e = 0;
e < n_element;
e++)
2458 for (
unsigned j = 0; j < n_internal; j++)
2474 const unsigned long n_node = this->
nnode();
2475 for (
unsigned long n = 0; n < n_node; n++)
2478 if ((this->
Node_pt[n]->does_pointer_correspond_to_value(parameter_pt)) ||
2479 (this->
Node_pt[n]->does_pointer_correspond_to_position_data(
2487 const unsigned long n_element = this->
nelement();
2488 for (
unsigned long e = 0;
e < n_element;
e++)
2497 for (
unsigned j = 0; j < n_internal; j++)
2517 const bool& preserve_existing_data)
2520 const unsigned long n_node = this->
nnode();
2521 for (
unsigned long n = 0; n < n_node; n++)
2524 this->
Node_pt[n]->set_time_stepper(time_stepper_pt,
2525 preserve_existing_data);
2526 this->
Node_pt[n]->set_position_time_stepper(time_stepper_pt,
2527 preserve_existing_data);
2536 TimeStepper*
const& time_stepper_pt,
const bool& preserve_existing_data)
2539 const unsigned long n_element = this->
nelement();
2540 for (
unsigned long e = 0;
e < n_element;
e++)
2543 const unsigned n_internal = this->
element_pt(
e)->ninternal_data();
2546 for (
unsigned j = 0; j < n_internal; j++)
2548 this->
element_pt(
e)->internal_data_pt(j)->set_time_stepper(
2549 time_stepper_pt, preserve_existing_data);
2568 for (
unsigned e = 0, ne =
nelement();
e < ne;
e++)
2605 std::list<std::pair<unsigned long, int>>
2606 list_of_elements_and_local_node_numbers;
2609 unsigned long n_element = this->
nelement();
2610 for (
unsigned long e = 0;
e < n_element;
e++)
2621 if (node_number != -1)
2623 list_of_elements_and_local_node_numbers.insert(
2624 list_of_elements_and_local_node_numbers.end(),
2625 std::make_pair(
e, node_number));
2633 if (list_of_elements_and_local_node_numbers.empty())
2635 std::ostringstream error_stream;
2636 error_stream <<
"Node " <<
node_pt
2637 <<
" is not contained in any elements in the Mesh."
2639 <<
"How was it created then?" << std::endl;
2642 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2653 std::list<std::pair<unsigned long, int>>::iterator list_it =
2654 list_of_elements_and_local_node_numbers.begin();
2667 if (solid_node_pt != 0)
2673 new_node_pt->
copy(old_node_pt);
2679 list_it != list_of_elements_and_local_node_numbers.end();
2705 #ifdef OOMPH_HAS_MPI
2721 double t_start = 0.0;
2728 int n_proc =
Comm_pt->nproc();
2729 int my_rank =
Comm_pt->my_rank();
2734 for (
int d = 0; d < n_proc; d++)
2738 std::map<Node*, bool> node_shared;
2747 unsigned nhalo_elem = halo_elem_pt.size();
2749 for (
unsigned e = 0;
e < nhalo_elem;
e++)
2756 unsigned nnod = el_pt->
nnode();
2757 for (
unsigned j = 0; j < nnod; j++)
2762 if (!node_shared[nod_pt])
2765 node_shared[nod_pt] =
true;
2775 unsigned nhaloed_elem = haloed_elem_pt.size();
2777 for (
unsigned e = 0;
e < nhaloed_elem;
e++)
2785 unsigned nnod = el_pt->
nnode();
2786 for (
unsigned j = 0; j < nnod; j++)
2791 if (!node_shared[nod_pt])
2794 node_shared[nod_pt] =
true;
2809 unsigned nhaloed_elem = haloed_elem_pt.size();
2811 for (
unsigned e = 0;
e < nhaloed_elem;
e++)
2819 unsigned nnod = el_pt->
nnode();
2820 for (
unsigned j = 0; j < nnod; j++)
2825 if (!node_shared[nod_pt])
2828 node_shared[nod_pt] =
true;
2837 unsigned nhalo_elem = halo_elem_pt.size();
2839 for (
unsigned e = 0;
e < nhalo_elem;
e++)
2846 unsigned nnod = el_pt->
nnode();
2847 for (
unsigned j = 0; j < nnod; j++)
2852 if (!node_shared[nod_pt])
2855 node_shared[nod_pt] =
true;
2871 oomph_info <<
"Time for identification of shared nodes: "
2872 << t_end - t_start << std::endl;
2892 double t_start = 0.0;
2898 double tt_start = 0.0;
2899 double tt_end = 0.0;
2906 int n_proc =
Comm_pt->nproc();
2907 int my_rank =
Comm_pt->my_rank();
2916 "Processor has shared nodes with itself! Something's gone wrong!",
2917 OOMPH_CURRENT_FUNCTION,
2918 OOMPH_EXCEPTION_LOCATION);
2925 std::map<Node*, std::set<int>> shared_domain_set;
2928 std::map<Node*, unsigned> global_haloed_node_number_plus_one;
2929 unsigned global_haloed_count = 0;
2932 for (
int d = 0; d < n_proc; d++)
2939 for (
unsigned j = 0; j < nnod_haloed; j++)
2942 shared_domain_set[nod_pt].insert(d);
2943 if (global_haloed_node_number_plus_one[nod_pt] == 0)
2945 global_haloed_node_number_plus_one[nod_pt] =
2946 global_haloed_count + 1;
2947 global_haloed_count++;
2956 oomph_info <<
"Time for initial classification in "
2957 "Mesh::synchronise_shared_nodes(): "
2958 << tt_end - tt_start << std::endl;
2980 for (
int domain = 0; domain < n_proc; domain++)
2983 send_displacement[domain] = send_data.size();
2988 if (domain != my_rank)
2995 for (
unsigned j = 0; j < nnod_haloed; j++)
3000 send_data.push_back(global_haloed_node_number_plus_one[nod_pt] - 1);
3003 std::set<int> tmp_shared_domain_set = shared_domain_set[nod_pt];
3006 unsigned n_shared_domains = tmp_shared_domain_set.size();
3007 send_data.push_back(n_shared_domains);
3010 for (std::set<int>::iterator it = tmp_shared_domain_set.begin();
3011 it != tmp_shared_domain_set.end();
3014 send_data.push_back(*it);
3020 send_n[domain] = send_data.size() - send_displacement[domain];
3029 &send_n[0], 1, MPI_INT, &receive_n[0], 1, MPI_INT,
Comm_pt->mpi_comm());
3035 int receive_data_count = 0;
3036 for (
int rank = 0; rank < n_proc; ++rank)
3039 receive_displacement[rank] = receive_data_count;
3040 receive_data_count += receive_n[rank];
3045 if (receive_data_count == 0)
3047 ++receive_data_count;
3053 if (send_data.size() == 0)
3055 send_data.resize(1);
3059 MPI_Alltoallv(&send_data[0],
3061 &send_displacement[0],
3065 &receive_displacement[0],
3073 oomph_info <<
"Time for alltoall in Mesh::synchronise_shared_nodes(): "
3074 << tt_end - tt_start << std::endl;
3081 oomph_info <<
"Starting vector to set conversion in "
3082 <<
"Mesh::synchronise_shared_nodes() for a total of "
3091 for (
int d = 0; d < n_proc; d++)
3094 for (
unsigned i = 0;
i < n_vector;
i++)
3105 <<
"Time for vector to set in Mesh::synchronise_shared_nodes(): "
3106 << tt_end - tt_start << std::endl;
3112 for (
int send_rank = 0; send_rank < n_proc; send_rank++)
3116 if ((send_rank != my_rank) && (receive_n[send_rank] != 0))
3119 unsigned count = receive_displacement[send_rank];
3129 std::map<unsigned, std::pair<Node*, std::set<unsigned>>> domains_map;
3132 unsigned nnod_halo = this->
nhalo_node(send_rank);
3133 for (
unsigned j = 0; j < nnod_halo; j++)
3138 unsigned haloed_node_id_on_send_proc = receive_data[count++];
3141 unsigned n_shared_domains = receive_data[count++];
3144 std::set<unsigned> domain_set;
3147 for (
unsigned i = 0;
i < n_shared_domains;
i++)
3149 int shared_domain = receive_data[count++];
3152 domain_set.insert(shared_domain);
3156 domains_map[haloed_node_id_on_send_proc] =
3157 std::make_pair(nod_pt, domain_set);
3164 int previous_one = -1;
3166 for (std::map<
unsigned, std::pair<
Node*, std::set<unsigned>>>::iterator
3167 it = domains_map.begin();
3168 it != domains_map.end();
3174 if (
int((*it).first) < previous_one)
3176 std::ostringstream error_stream;
3177 error_stream <<
"Map did not store entries in order of key\n "
3178 <<
"Current key: " << (*it).first
3179 <<
"; previous one: " << previous_one <<
"\n"
3180 <<
"Need to rewrite this...\n";
3182 OOMPH_CURRENT_FUNCTION,
3183 OOMPH_EXCEPTION_LOCATION);
3185 previous_one = (*it).first;
3189 Node* nod_pt = (*it).second.first;
3192 std::set<unsigned> domain_set((*it).second.second);
3195 for (std::set<unsigned>::iterator itt = domain_set.begin();
3196 itt != domain_set.end();
3199 int shared_domain = (*itt);
3202 if (shared_domain != my_rank)
3206 std::set<Node*>::iterator ittt =
3207 shared_node_set[shared_domain].find(nod_pt);
3211 if (ittt == shared_node_set[shared_domain].end())
3217 shared_node_set[shared_domain].insert(nod_pt);
3234 "Processor has shared nodes with itself! Something's gone wrong!",
3235 OOMPH_CURRENT_FUNCTION,
3236 OOMPH_EXCEPTION_LOCATION);
3245 <<
"Time for final processing in Mesh::synchronise_shared_nodes(): "
3246 << tt_end - tt_start << std::endl;
3253 oomph_info <<
"Total time for Mesh::synchronise_shared_nodes(): "
3254 << t_end - t_start << std::endl;
3263 const bool& report_stats)
3268 double t_start = 0.0;
3274 double tt_start = 0.0;
3285 double tt_end = 0.0;
3289 oomph_info <<
"Time for Mesh::setup_shared_node_scheme() "
3290 <<
" Mesh::classify_halo_and_haloed_nodes(): "
3291 << tt_end - tt_start << std::endl;
3302 int n_proc =
Comm_pt->nproc();
3303 int my_rank =
Comm_pt->my_rank();
3308 std::map<Data*, std::set<unsigned>> processors_associated_with_data;
3309 std::map<Data*, unsigned> processor_in_charge;
3313 for (
int domain = 0; domain < n_proc; domain++)
3319 unsigned nelem = halo_elem_pt.size();
3320 for (
unsigned e = 0;
e < nelem;
e++)
3325 if (finite_el_pt != 0)
3328 unsigned nnod = finite_el_pt->
nnode();
3329 for (
unsigned j = 0; j < nnod; j++)
3333 processors_associated_with_data[nod_pt].insert(domain);
3337 if (solid_nod_pt != 0)
3339 processors_associated_with_data[solid_nod_pt
3352 for (
unsigned e = 0;
e < nelem;
e++)
3358 if ((finite_el_pt != 0) && (!finite_el_pt->
is_halo()))
3361 unsigned nnod = finite_el_pt->
nnode();
3362 for (
unsigned j = 0; j < nnod; j++)
3367 processors_associated_with_data[nod_pt].insert(my_rank);
3371 if (solid_nod_pt != 0)
3373 processors_associated_with_data[solid_nod_pt
3386 <<
"Time for setup loops in Mesh::classify_halo_and_haloed_nodes: "
3387 << tt_end - tt_start << std::endl;
3403 for (
int d = 0; d < n_proc; d++)
3416 unsigned count_data = 0;
3419 unsigned n_haloed_elem = haloed_elem_pt.size();
3420 for (
unsigned e = 0;
e < n_haloed_elem;
e++)
3425 if (haloed_el_pt != 0)
3428 unsigned n_node = haloed_el_pt->
nnode();
3429 for (
unsigned j = 0; j < n_node; j++)
3434 unsigned n_assoc = processors_associated_with_data[nod_pt].size();
3437 processors_associated_with_data_on_other_proc.push_back(n_assoc);
3441 std::set<unsigned> procs_set =
3442 processors_associated_with_data[nod_pt];
3443 for (std::set<unsigned>::iterator it = procs_set.begin();
3444 it != procs_set.end();
3447 processors_associated_with_data_on_other_proc.push_back(*it);
3456 MPI_Send(&count_data, 1, MPI_UNSIGNED, d, 0,
Comm_pt->mpi_comm());
3457 if (count_data != 0)
3459 MPI_Send(&processors_associated_with_data_on_other_proc[0],
3470 for (
int dd = 0; dd < n_proc; dd++)
3476 unsigned n_halo_elem = halo_elem_pt.size();
3477 unsigned count_data = 0;
3478 MPI_Recv(&count_data,
3485 if (count_data != 0)
3487 processors_associated_with_data_on_other_proc.resize(count_data);
3488 MPI_Recv(&processors_associated_with_data_on_other_proc[0],
3498 for (
unsigned e = 0;
e < n_halo_elem;
e++)
3502 if (halo_el_pt != 0)
3504 unsigned n_node = halo_el_pt->
nnode();
3505 for (
unsigned j = 0; j < n_node; j++)
3512 processors_associated_with_data_on_other_proc[count_data];
3515 for (
unsigned i_assoc = 0; i_assoc < n_assoc; i_assoc++)
3518 unsigned sent_domain =
3519 processors_associated_with_data_on_other_proc
3524 processors_associated_with_data[nod_pt].insert(
3530 if (solid_nod_pt != 0)
3532 processors_associated_with_data
3534 .insert(sent_domain);
3550 <<
"Time for pt2pt send/recv in Mesh::classify_halo_and_haloed_nodes: "
3551 << tt_end - tt_start << std::endl;
3561 unsigned nnod = this->
nnode();
3562 for (
unsigned j = 0; j < nnod; j++)
3572 if (solid_nod_pt != 0)
3578 unsigned proc_max = 0;
3579 std::set<unsigned> procs_set = processors_associated_with_data[nod_pt];
3580 for (std::set<unsigned>::iterator it = procs_set.begin();
3581 it != procs_set.end();
3584 if (*it > proc_max) proc_max = *it;
3586 processor_in_charge[nod_pt] = proc_max;
3589 if (solid_nod_pt != 0)
3592 unsigned proc_max_solid = 0;
3593 std::set<unsigned> procs_set_solid =
3595 for (std::set<unsigned>::iterator it = procs_set_solid.begin();
3596 it != procs_set_solid.end();
3599 if (*it > proc_max_solid) proc_max_solid = *it;
3612 std::map<Node*, bool> done;
3615 for (
int domain = 0; domain < n_proc; domain++)
3621 unsigned nelem = halo_elem_pt.size();
3623 for (
unsigned e = 0;
e < nelem;
e++)
3630 if (finite_el_pt != 0)
3633 unsigned nnod = finite_el_pt->
nnode();
3634 for (
unsigned j = 0; j < nnod; j++)
3642 int proc_in_charge = processor_in_charge[nod_pt];
3644 if (proc_in_charge != my_rank)
3649 if (proc_in_charge ==
int(domain))
3661 if (solid_nod_pt != 0)
3668 done[nod_pt] =
true;
3678 for (
unsigned iintern = 0; iintern < nintern_data; iintern++)
3690 for (
int domain = 0; domain < n_proc; domain++)
3693 std::map<Node*, bool> node_done;
3700 unsigned nelem = haloed_elem_pt.size();
3702 for (
unsigned e = 0;
e < nelem;
e++)
3709 if (finite_el_pt != 0)
3712 unsigned nnod = finite_el_pt->
nnode();
3713 for (
unsigned j = 0; j < nnod; j++)
3718 if (!node_done[nod_pt])
3721 int proc_in_charge = processor_in_charge[nod_pt];
3723 if (proc_in_charge == my_rank)
3729 node_done[nod_pt] =
true;
3742 <<
"Time for first classific in Mesh::classify_halo_and_haloed_nodes: "
3743 << tt_end - tt_start << std::endl;
3770 unsigned n_overlooked_halo = 0;
3790 for (
int domain = 0; domain < n_proc; domain++)
3793 send_displacement[domain] = send_data.size();
3797 if (domain != my_rank)
3800 for (
unsigned j = 0; j < nnod; j++)
3805 int proc_in_charge = processor_in_charge[nod_pt];
3806 if ((proc_in_charge != my_rank) && !(nod_pt->
is_halo()))
3813 n_overlooked_halo++;
3814 over_looked_halo_node_pt[proc_in_charge].push_back(nod_pt);
3822 if (solid_nod_pt != 0)
3830 send_data.push_back(j);
3831 send_data.push_back(proc_in_charge);
3837 send_data.push_back(-1);
3840 send_n[domain] = send_data.size() - send_displacement[domain];
3846 unsigned global_max_n_overlooked_halo = 0;
3847 MPI_Allreduce(&n_overlooked_halo,
3848 &global_max_n_overlooked_halo,
3855 oomph_info <<
"Global max number of overlooked haloes: "
3856 << global_max_n_overlooked_halo << std::endl;
3861 oomph_info <<
"Time for setup 1st alltoalls in "
3862 "Mesh::classify_halo_and_haloed_nodes: "
3863 << tt_end - tt_start << std::endl;
3871 if (global_max_n_overlooked_halo > 0)
3878 &send_n[0], 1, MPI_INT, &receive_n[0], 1, MPI_INT,
Comm_pt->mpi_comm());
3885 <<
"Time for 1st alltoall in Mesh::classify_halo_and_haloed_nodes: "
3886 << tt_end - tt_start << std::endl;
3897 int receive_data_count = 0;
3898 for (
int rank = 0; rank < n_proc; ++rank)
3901 receive_displacement[rank] = receive_data_count;
3902 receive_data_count += receive_n[rank];
3907 if (receive_data_count == 0)
3909 ++receive_data_count;
3915 if (send_data.size() == 0)
3917 send_data.resize(1);
3921 MPI_Alltoallv(&send_data[0],
3923 &send_displacement[0],
3927 &receive_displacement[0],
3936 <<
"Time for 2nd alltoall in Mesh::classify_halo_and_haloed_nodes: "
3937 << tt_end - tt_start << std::endl;
3949 for (
int send_rank = 0; send_rank < n_proc; send_rank++)
3953 if ((send_rank != my_rank) && (receive_n[send_rank] != 0))
3956 unsigned count = receive_displacement[send_rank];
3962 int next_one = receive_data[count++];
3972 unsigned j = unsigned(next_one);
3975 unsigned proc_in_charge = unsigned(receive_data[count++]);
3991 for (
unsigned jj = 0; jj < nnod; jj++)
3999 send_data_for_proc_in_charge[proc_in_charge].push_back(jj);
4002 send_data_for_proc_in_charge[proc_in_charge].push_back(
4010 std::ostringstream error_stream;
4012 <<
"Failed to find node that is shared node " << j
4013 <<
" (with processor " << send_rank
4014 <<
") \n in shared node lookup scheme with processor "
4015 << proc_in_charge <<
" which is in charge.\n";
4017 OOMPH_CURRENT_FUNCTION,
4018 OOMPH_EXCEPTION_LOCATION);
4029 oomph_info <<
"Time for 1st setup 3rd alltoall in "
4030 "Mesh::classify_halo_and_haloed_nodes: "
4031 << tt_end - tt_start << std::endl;
4048 for (
int domain = 0; domain < n_proc; domain++)
4051 all_send_displacement[domain] = all_send_data.size();
4055 if (domain != my_rank)
4057 unsigned n = send_data_for_proc_in_charge[domain].size();
4058 for (
unsigned j = 0; j < n; j++)
4060 all_send_data.push_back(send_data_for_proc_in_charge[domain][j]);
4065 all_send_data.push_back(-1);
4068 all_send_n[domain] =
4069 all_send_data.size() - all_send_displacement[domain];
4076 oomph_info <<
"Time for 2nd setup 3rd alltoall in "
4077 "Mesh::classify_halo_and_haloed_nodes: "
4078 << tt_end - tt_start << std::endl;
4089 MPI_Alltoall(&all_send_n[0],
4102 <<
"Time for 3rd alltoall in Mesh::classify_halo_and_haloed_nodes: "
4103 << tt_end - tt_start << std::endl;
4113 int all_receive_data_count = 0;
4115 for (
int rank = 0; rank < n_proc; ++rank)
4118 all_receive_displacement[rank] = all_receive_data_count;
4119 all_receive_data_count += all_receive_n[rank];
4124 if (all_receive_data_count == 0)
4126 ++all_receive_data_count;
4128 Vector<int> all_receive_data(all_receive_data_count);
4132 if (all_send_data.size() == 0)
4134 all_send_data.resize(1);
4138 MPI_Alltoallv(&all_send_data[0],
4140 &all_send_displacement[0],
4142 &all_receive_data[0],
4144 &all_receive_displacement[0],
4153 <<
"Time for 4th alltoall in Mesh::classify_halo_and_haloed_nodes: "
4154 << tt_end - tt_start << std::endl;
4163 for (
int send_rank = 0; send_rank < n_proc; send_rank++)
4167 if ((send_rank != my_rank) && (all_receive_n[send_rank] != 0))
4170 unsigned count = all_receive_displacement[send_rank];
4176 int next_one = all_receive_data[count++];
4186 unsigned j = unsigned(next_one);
4192 unsigned proc_with_overlooked_halo = all_receive_data[count++];
4199 over_looked_haloed_node_pt[proc_with_overlooked_halo].push_back(
4209 oomph_info <<
"Time for postproc 4th alltoall in "
4210 "Mesh::classify_halo_and_haloed_nodes: "
4211 << tt_end - tt_start << std::endl;
4220 for (
int d = 0; d < n_proc; d++)
4227 unsigned nnod = over_looked_halo_node_pt[d].size();
4228 for (
unsigned j = 0; j < nnod; j++)
4232 nnod = over_looked_haloed_node_pt[d].size();
4233 for (
unsigned j = 0; j < nnod; j++)
4238 else if (d > my_rank)
4240 unsigned nnod = over_looked_haloed_node_pt[d].size();
4241 for (
unsigned j = 0; j < nnod; j++)
4245 nnod = over_looked_halo_node_pt[d].size();
4246 for (
unsigned j = 0; j < nnod; j++)
4258 oomph_info <<
"BEFORE SYNCHRONISE SHARED NODES Processor " << my_rank
4259 <<
" holds " << this->
nnode() <<
" nodes of which "
4260 << this->
nhalo_node() <<
" are halo nodes \n while "
4262 << this->
nshared_node() <<
" are shared nodes." << std::endl;
4266 for (
int iproc = 0; iproc < n_proc; iproc++)
4273 oomph_info <<
"With process " << iproc <<
", there are "
4274 << this->
nhalo_node(iproc) <<
" halo nodes, and "
4279 << halo_elem_pt.size() <<
" halo elements and "
4280 << haloed_elem_pt.size() <<
" haloed elements\n";
4331 "Processor has haloed nodes with itself! Something's gone wrong!",
4332 OOMPH_CURRENT_FUNCTION,
4333 OOMPH_EXCEPTION_LOCATION);
4339 "Processor has halo nodes with itself! Something's gone wrong!",
4340 OOMPH_CURRENT_FUNCTION,
4341 OOMPH_EXCEPTION_LOCATION);
4346 throw OomphLibError(
"Processor has root haloed elements with itself! "
4347 "Something's gone wrong!",
4348 OOMPH_CURRENT_FUNCTION,
4349 OOMPH_EXCEPTION_LOCATION);
4355 "Processor has root halo elements with itself! Something's gone wrong!",
4356 OOMPH_CURRENT_FUNCTION,
4357 OOMPH_EXCEPTION_LOCATION);
4369 <<
" are shared nodes." << std::endl;
4373 for (
int iproc = 0; iproc < n_proc; iproc++)
4379 oomph_info <<
"With process " << iproc <<
", there are "
4380 << this->
nhalo_node(iproc) <<
" halo nodes, and "
4383 << this->
nshared_node(iproc) <<
" shared nodes" << std::endl
4384 << halo_elem_pt.size() <<
" halo elements and "
4385 << haloed_elem_pt.size() <<
" haloed elements\n";
4406 oomph_info <<
"Total time for Mesh::classify_halo_and_halo_nodes(): "
4407 << t_end - t_start << std::endl;
4428 double t_start = 0.0;
4440 int n_proc =
Comm_pt->nproc();
4441 int my_rank =
Comm_pt->my_rank();
4444 for (
int d = 0; d < n_proc; d++)
4451 for (
int dd = 0; dd < n_proc; dd++)
4460 if ((nnod_halo + nnod_ext_halo) != 0)
4467 &tmp[0], 3, MPI_INT, dd, 0,
Comm_pt->mpi_comm(), &status);
4471 int nnod_haloed = tmp[0];
4472 if (nnod_haloed != nnod_halo)
4474 std::ostringstream error_message;
4475 error_message <<
"Clash in numbers of halo and haloed nodes "
4477 error_message <<
" between procs " << dd <<
" and " << d
4478 <<
": " << nnod_haloed <<
" " << nnod_halo
4481 OOMPH_CURRENT_FUNCTION,
4482 OOMPH_EXCEPTION_LOCATION);
4486 int nnod_ext_haloed = tmp[1];
4487 if (nnod_ext_haloed != nnod_ext_halo)
4489 std::ostringstream error_message;
4491 <<
"Clash in numbers of external halo and haloed nodes "
4493 error_message <<
" between procs " << dd <<
" and " << d
4494 <<
": " << nnod_ext_haloed <<
" "
4495 << nnod_ext_halo << std::endl;
4497 OOMPH_CURRENT_FUNCTION,
4498 OOMPH_EXCEPTION_LOCATION);
4503 unsigned n_rec = tmp[2];
4507 MPI_Recv(&unsigned_rec_data[0],
4519 for (
unsigned loop = 0; loop < 2; loop++)
4521 unsigned hi_nod = nnod_halo;
4524 hi_nod = nnod_ext_halo;
4526 for (
int j = 0; j < int(hi_nod); j++)
4540 unsigned nval_local = nod_pt->
nvalue();
4543 unsigned nval_other = unsigned_rec_data[count++];
4545 if (nval_local != nval_other)
4547 nod_pt->
resize(nval_other);
4551 unsigned nentry = unsigned_rec_data[count++];
4561 "Failed to cast node to boundary node even though "
4562 "we've received data for boundary node",
4563 OOMPH_CURRENT_FUNCTION,
4564 OOMPH_EXCEPTION_LOCATION);
4569 bool already_existed =
true;
4572 ->index_of_first_value_assigned_by_face_element_pt() ==
4577 new std::map<unsigned, unsigned>;
4578 already_existed =
false;
4583 std::map<unsigned, unsigned>* map_pt =
4588 for (
unsigned i = 0;
i < nentry;
i++)
4591 unsigned first_received = unsigned_rec_data[count++];
4592 unsigned second_received = unsigned_rec_data[count++];
4595 if (already_existed)
4598 if ((*map_pt)[first_received] != second_received)
4600 std::ostringstream error_message;
4601 error_message <<
"Existing map entry for map entry "
4602 <<
i <<
" for node located at ";
4603 unsigned n = nod_pt->
ndim();
4604 for (
unsigned ii = 0; ii < n; ii++)
4606 error_message << nod_pt->
position(ii) <<
" ";
4609 <<
"Key: " << first_received <<
" "
4610 <<
"Local value: " << (*map_pt)[first_received]
4612 <<
"Received value: " << second_received
4615 OOMPH_CURRENT_FUNCTION,
4616 OOMPH_EXCEPTION_LOCATION);
4623 (*map_pt)[first_received] = second_received;
4642 tmp[0] = nnod_haloed;
4644 tmp[1] = nnod_ext_haloed;
4645 if ((nnod_haloed + nnod_ext_haloed) != 0)
4651 for (
unsigned loop = 0; loop < 2; loop++)
4653 unsigned hi_nod = nnod_haloed;
4656 hi_nod = nnod_ext_haloed;
4658 for (
int j = 0; j < int(hi_nod); j++)
4672 unsigned_send_data.push_back(nod_pt->
nvalue());
4683 unsigned_send_data.push_back(0);
4688 std::map<unsigned, unsigned>* map_pt =
4695 unsigned_send_data.push_back(0);
4701 unsigned_send_data.push_back(map_pt->size());
4704 for (std::map<unsigned, unsigned>::iterator p =
4709 unsigned_send_data.push_back((*p).first);
4710 unsigned_send_data.push_back((*p).second);
4718 int n_send = unsigned_send_data.size();
4722 MPI_Send(&tmp[0], 3, MPI_INT, d, 0,
Comm_pt->mpi_comm());
4725 MPI_Send(&unsigned_send_data[0],
4739 oomph_info <<
"Total time for Mesh::resize_halo_nodes(): "
4740 << t_end - t_start << std::endl;
4757 unsigned n_node = (it->second).size();
4759 for (
unsigned n = 0; n < n_node; n++)
4763 (it->second)[n]->add_value_pt_to_map(map_of_halo_data);
4775 unsigned n_element = (it->second).size();
4776 for (
unsigned e = 0;
e < n_element;
e++)
4791 unsigned nleaf = leaf_pt.size();
4792 for (
unsigned l = 0; l < nleaf; l++)
4794 leaf_pt[l]->object_pt()->add_internal_value_pt_to_map(
4812 unsigned n_node = (it->second).size();
4814 for (
unsigned n = 0; n < n_node; n++)
4818 (it->second)[n]->add_value_pt_to_map(map_of_halo_data);
4830 unsigned n_element = (it->second).size();
4831 for (
unsigned e = 0;
e < n_element;
e++)
4833 (it->second)[
e]->add_internal_value_pt_to_map(map_of_halo_data);
4846 unsigned& max_number,
4847 unsigned& min_number)
4850 int n_proc =
Comm_pt->nproc();
4851 int my_rank =
Comm_pt->my_rank();
4863 MPI_Gather(&nhalo_node_local,
4875 int min = 1000000000;
4879 for (
int i = 0;
i < n_proc;
i++)
4881 av_number += double(nhalo_nodes[
i]);
4882 if (
int(nhalo_nodes[
i]) > max) max = nhalo_nodes[
i];
4883 if (
int(nhalo_nodes[
i]) < min) min = nhalo_nodes[
i];
4885 av_number /= double(n_proc);
4889 MPI_Bcast(&max, 1, MPI_INT, 0,
Comm_pt->mpi_comm());
4890 MPI_Bcast(&min, 1, MPI_INT, 0,
Comm_pt->mpi_comm());
4891 MPI_Bcast(&av_number, 1, MPI_DOUBLE, 0,
Comm_pt->mpi_comm());
4905 unsigned& max_number,
4906 unsigned& min_number)
4909 int n_proc =
Comm_pt->nproc();
4910 int my_rank =
Comm_pt->my_rank();
4922 MPI_Gather(&nhaloed_node_local,
4934 int min = 1000000000;
4938 for (
int i = 0;
i < n_proc;
i++)
4940 av_number += double(nhaloed_nodes[
i]);
4941 if (
int(nhaloed_nodes[
i]) > max) max = nhaloed_nodes[
i];
4942 if (
int(nhaloed_nodes[
i]) < min) min = nhaloed_nodes[
i];
4944 av_number /= double(n_proc);
4948 MPI_Bcast(&max, 1, MPI_INT, 0,
Comm_pt->mpi_comm());
4949 MPI_Bcast(&min, 1, MPI_INT, 0,
Comm_pt->mpi_comm());
4950 MPI_Bcast(&av_number, 1, MPI_DOUBLE, 0,
Comm_pt->mpi_comm());
4963 const bool& report_stats,
4964 const bool& overrule_keep_as_halo_element_status)
4970 int n_proc = comm_pt->nproc();
4971 int my_rank = comm_pt->my_rank();
4975 unsigned nnod = this->
nnode();
4977 std::ostringstream filename;
4987 for (
int d = 0; d < n_proc; d++)
4993 filename << doc_info.
directory() <<
"/domain" << d <<
"-"
4994 << doc_info.
number() <<
".dat";
4995 domain_file[d] =
new std::ofstream(filename.str().c_str());
4999 for (
unsigned e = 0;
e < nelem;
e++)
5007 f_el_pt->
output(*domain_file[element_domain[
e]], 5);
5011 for (
int d = 0; d < n_proc; d++)
5013 domain_file[d]->close();
5014 delete domain_file[d];
5028 std::map<Data*, std::set<unsigned>> processors_associated_with_data;
5029 std::map<Data*, unsigned> processor_in_charge;
5032 for (
unsigned j = 0; j < nnod; j++)
5035 processor_in_charge[nod_pt] = 0;
5039 for (
unsigned e = 0;
e < nelem;
e++)
5046 unsigned el_domain = element_domain[
e];
5049 unsigned nnod = el_pt->
nnode();
5050 for (
unsigned j = 0; j < nnod; j++)
5055 if (el_domain > processor_in_charge[nod_pt])
5057 processor_in_charge[nod_pt] = el_domain;
5059 processors_associated_with_data[nod_pt].insert(el_domain);
5072 for (
int d = 0; d < n_proc; d++)
5078 filename << doc_info.
directory() <<
"/node" << d <<
"-"
5079 << doc_info.
number() <<
".dat";
5080 node_file[d] =
new std::ofstream(filename.str().c_str());
5084 for (
unsigned j = 0; j < nnod; j++)
5087 const unsigned n_dim = nod_pt->
ndim();
5088 for (
unsigned i = 0;
i < n_dim;
i++)
5090 *node_file[processor_in_charge[nod_pt]] << nod_pt->
x(
i) <<
" ";
5092 *node_file[processor_in_charge[nod_pt]] <<
"\n";
5094 for (
int d = 0; d < n_proc; d++)
5096 node_file[d]->close();
5097 delete node_file[d];
5106 for (
unsigned j = 0; j < nnod; j++)
5117 for (
unsigned e = 0;
e < nelem;
e++)
5127 const unsigned tmp_nboundary = this->
nboundary();
5138 bool is_a_triangle_mesh_base_mesh =
false;
5139 if (triangle_mesh_pt != 0)
5143 is_a_triangle_mesh_base_mesh =
true;
5146 const unsigned n_regions = triangle_mesh_pt->
nregion();
5149 for (
unsigned ib = 0; ib < tmp_nboundary; ib++)
5155 ntmp_boundary_elements_in_region[ib].resize(n_regions);
5158 for (
unsigned rr = 0; rr < n_regions; rr++)
5161 const unsigned region_id =
5167 ntmp_boundary_elements_in_region[ib][rr] =
5174 for (
unsigned e = 0;
e < nelem;
e++)
5190 std::vector<bool> element_retained(nelem,
false);
5200 for (
int p = 0; p < n_proc; p++)
5202 root_halo_element[p].push_back(-1);
5207 unsigned number_of_retained_elements = 0;
5210 nelem = backed_up_el_pt.size();
5211 for (
unsigned e = 0;
e < nelem;
e++)
5215 unsigned el_domain = element_domain[
e];
5218 if (el_domain ==
unsigned(my_rank))
5221 element_retained[
e] =
true;
5222 number_of_retained_elements++;
5233 if (!overrule_keep_as_halo_element_status)
5237 if (!element_retained[
e])
5239 root_halo_element[el_domain].push_back(
e);
5240 element_retained[
e] =
true;
5241 number_of_retained_elements++;
5250 if (finite_el_pt != 0)
5252 unsigned n_node = finite_el_pt->
nnode();
5253 for (
unsigned n = 0; n < n_node; n++)
5258 unsigned keep_it =
false;
5259 for (std::set<unsigned>::iterator it =
5260 processors_associated_with_data[nod_pt].begin();
5261 it != processors_associated_with_data[nod_pt].end();
5264 if (*it ==
unsigned(my_rank))
5277 if (!element_retained[
e])
5279 root_halo_element[el_domain].push_back(
e);
5280 element_retained[
e] =
true;
5281 number_of_retained_elements++;
5300 if (is_a_triangle_mesh_base_mesh)
5313 processors_associated_with_data,
5314 overrule_keep_as_halo_element_status);
5331 nelem = backed_up_el_pt.size();
5332 for (
unsigned e = 0;
e < nelem;
e++)
5335 if (element_retained[
e])
5354 if (!is_a_triangle_mesh_base_mesh)
5356 deleted_element_pt.push_back(el_pt);
5359 if (is_a_triangle_mesh_base_mesh)
5362 deleted_f_element_pt.push_back(backed_up_f_el_pt[
e]);
5376 std::map<unsigned, bool> done;
5378 for (
int d = 0; d < n_proc; d++)
5380 nelem = root_halo_element[d].size();
5381 for (
unsigned e = 0;
e < nelem;
e++)
5383 int number = root_halo_element[d][
e];
5389 std::ostringstream error_message;
5390 error_message <<
"Have already added element " << number
5391 <<
" as root halo element\n"
5394 OOMPH_CURRENT_FUNCTION,
5395 OOMPH_EXCEPTION_LOCATION);
5397 done[number] =
true;
5415 for (
int p = 0; p < n_proc; p++)
5417 nhalo[p] = root_halo_element[p].size();
5422 for (
int p = 0; p < n_proc; p++)
5428 &nhalo[p], 1, MPI_INT, &nhaloed[0], 1, MPI_INT, p, comm_pt->mpi_comm());
5436 unsigned total_number_of_root_haloed_elements = 0;
5437 for (
int i_proc = 0; i_proc < n_proc; i_proc++)
5439 total_number_of_root_haloed_elements += nhaloed[i_proc];
5442 start_index[i_proc] =
5443 total_number_of_root_haloed_elements - nhaloed[i_proc];
5454 Vector<int> all_root_haloed_element(total_number_of_root_haloed_elements,
5458 for (
int p = 0; p < n_proc; p++)
5463 MPI_Gatherv(&root_halo_element[p][0],
5467 &all_root_haloed_element[0],
5479 comm_pt->mpi_comm());
5488 for (
int d = 0; d < n_proc; d++)
5491 std::map<unsigned, bool> done;
5495 unsigned n = nhaloed[d];
5496 for (
unsigned e = 0;
e < n;
e++)
5498 int number = all_root_haloed_element[count];
5514 std::ostringstream error_message;
5515 error_message <<
"Have already added element " << number
5516 <<
" as root haloed element\n"
5519 OOMPH_CURRENT_FUNCTION,
5520 OOMPH_EXCEPTION_LOCATION);
5522 done[number] =
true;
5538 <<
" are root halo elements \n while "
5549 for (
unsigned e = 0;
e < nelem;
e++)
5558 unsigned nnod = f_el_pt->
nnode();
5559 for (
unsigned j = 0; j < nnod; j++)
5572 #ifdef OOMPH_HAS_TRIANGLE_LIB
5573 if (is_a_triangle_mesh_base_mesh)
5576 ntmp_boundary_elements,
5577 ntmp_boundary_elements_in_region,
5578 deleted_f_element_pt);
5587 #ifdef OOMPH_HAS_TRIANGLE_LIB
5596 if (ref_mesh_pt != 0)
5624 const bool& report_stats)
5629 #ifdef OOMPH_HAS_MPI
5637 if (ref_mesh_pt != 0)
5640 int n_proc =
Comm_pt->nproc();
5641 int my_rank =
Comm_pt->my_rank();
5647 <<
"\n----------------------------------------------------\n";
5648 oomph_info <<
"Before pruning: Processor " << my_rank <<
" holds "
5649 << this->
nelement() <<
" elements of which "
5651 <<
" are root halo elements \n while "
5653 <<
" are root haloed elements" << std::endl;
5656 oomph_info <<
"Before pruning: Processor " << my_rank <<
" holds "
5660 <<
" are shared nodes." << std::endl;
5664 for (
int iproc = 0; iproc < n_proc; iproc++)
5666 oomph_info <<
"Before pruning: With process " << iproc
5668 <<
" halo nodes, and " << std::endl
5674 <<
"----------------------------------------------------\n\n";
5678 double t_start = 0.0;
5687 unsigned nnod = this->
nnode();
5688 for (
unsigned j = 0; j < nnod; j++)
5697 std::map<RefineableElement*, bool> keep_element;
5698 for (
unsigned e = 0;
e < nelem;
e++)
5706 unsigned min_ref = 0;
5707 unsigned max_ref = 0;
5710 unsigned local_min_ref = 0;
5711 unsigned local_max_ref = 0;
5717 int int_local_min_ref = local_min_ref;
5718 int int_local_max_ref = local_max_ref;
5722 int_local_min_ref = INT_MAX;
5723 int_local_max_ref = INT_MIN;
5726 int int_min_ref = 0;
5727 int int_max_ref = 0;
5729 MPI_Allreduce(&int_local_min_ref,
5735 MPI_Allreduce(&int_local_max_ref,
5742 min_ref = unsigned(int_min_ref);
5743 max_ref = unsigned(int_max_ref);
5751 int local_n_ref = current_refined.size();
5756 local_n_ref = INT_MIN;
5762 &local_n_ref, &int_n_ref, 1, MPI_INT, MPI_MAX,
Comm_pt->mpi_comm());
5763 unsigned n_ref = int(int_n_ref);
5771 <<
"Time for establishing refinement levels in "
5772 <<
" Mesh::prune_halo_elements_and_nodes() [includes comms]: "
5773 << t_end - t_start << std::endl;
5784 unsigned base_level = n_ref - (max_ref - min_ref);
5793 base_level_elements_pt);
5796 unsigned n_base_el = base_level_elements_pt.size();
5797 for (
unsigned e = 0;
e < n_base_el;
e++)
5809 keep_element[ref_el_pt] =
true;
5812 unsigned nnod = ref_el_pt->
nnode();
5813 for (
unsigned j = 0; j < nnod; j++)
5826 unsigned n_unif = 0;
5827 unsigned n_unif_local =
5830 &n_unif_local, &n_unif, 1, MPI_UNSIGNED, MPI_MAX,
Comm_pt->mpi_comm());
5839 <<
"Time for synchronising refinement levels in "
5840 <<
" Mesh::prune_halo_elements_and_nodes() [includes comms]: "
5841 << t_end - t_start << std::endl;
5856 std::map<unsigned, Vector<FiniteElement*>> tmp_root_halo_element_pt;
5859 std::map<unsigned, Vector<FiniteElement*>> tmp_root_haloed_element_pt;
5862 std::map<unsigned, std::map<FiniteElement*, bool>>
5863 tmp_root_halo_element_is_retained;
5866 std::map<FiniteElement*, bool> halo_element_is_retained;
5868 for (
int domain = 0; domain < n_proc; domain++)
5874 unsigned nelem = halo_elem_pt.size();
5875 for (
unsigned e = 0;
e < nelem;
e++)
5886 if (halo_el_level == min_ref)
5896 el_pt = father_el_pt;
5903 unsigned nnod = el_pt->
nnode();
5904 for (
unsigned j = 0; j < nnod; j++)
5912 if (!tmp_root_halo_element_is_retained[domain][el_pt])
5914 tmp_root_halo_element_pt[domain].push_back(el_pt);
5915 tmp_root_halo_element_is_retained[domain][el_pt] =
true;
5917 keep_element[el_pt] =
true;
5918 halo_element_is_retained[el_pt] =
true;
5928 oomph_info <<
"Time for setup of retention pattern in "
5929 <<
" Mesh::prune_halo_elements_and_nodes(): "
5930 << t_end - t_start << std::endl;
5937 MPI_Barrier(
Comm_pt->mpi_comm());
5947 for (
int d = 0; d < n_proc; d++)
5950 for (
int dd = 0; dd < n_proc; dd++)
5964 unsigned nelem = haloed_elem_pt.size();
5971 MPI_Recv(&halo_kept[0],
5980 for (
unsigned e = 0;
e < nelem;
e++)
5992 if (haloed_el_level == min_ref)
6003 el_pt = father_el_pt;
6006 if (halo_kept[
e] == 1)
6010 bool already_root_haloed =
false;
6011 unsigned n_root_haloed =
6012 tmp_root_haloed_element_pt[dd].
size();
6013 for (
unsigned e_root = 0; e_root < n_root_haloed; e_root++)
6015 if (el_pt == tmp_root_haloed_element_pt[dd][e_root])
6017 already_root_haloed =
true;
6021 if (!already_root_haloed)
6023 tmp_root_haloed_element_pt[dd].push_back(el_pt);
6041 unsigned nelem = halo_elem_pt.size();
6043 for (
unsigned e = 0;
e < nelem;
e++)
6054 if (halo_el_level == min_ref)
6065 el_pt = father_el_pt;
6068 if (halo_element_is_retained[el_pt])
6080 &halo_kept[0], nelem, MPI_INT, d, 0,
Comm_pt->mpi_comm());
6091 oomph_info <<
"Time for pt2pt comms of retention pattern in "
6092 <<
" Mesh::prune_halo_elements_and_nodes(): "
6093 << t_end - t_start << std::endl;
6102 nnod = this->
nnode();
6104 for (
unsigned j = 0; j < nnod; j++)
6106 backed_up_nod_pt[j] = this->
node_pt(j);
6113 std::set<Tree*> trees_to_be_deleted_pt;
6116 nelem = backed_up_el_pt.size();
6117 for (
unsigned e = 0;
e < nelem;
e++)
6128 if (level == min_ref)
6138 el_pt = father_el_pt;
6143 if (keep_element[el_pt])
6163 tmp_all_tree_nodes_pt);
6167 unsigned n_tree = tmp_all_tree_nodes_pt.size();
6168 for (
unsigned j = 0; j < n_tree; j++)
6170 if (tmp_all_tree_nodes_pt[j]->object_pt() != 0)
6173 tmp_all_tree_nodes_pt[j]->object_pt()->refinement_level();
6176 if (!keep_element[tmp_all_tree_nodes_pt[j]->object_pt()])
6178 trees_to_be_deleted_pt.insert(tmp_all_tree_nodes_pt[j]);
6185 deleted_element_pt.push_back(ref_el_pt);
6193 for (std::set<Tree*>::iterator it = trees_to_be_deleted_pt.begin();
6194 it != trees_to_be_deleted_pt.end();
6197 Tree* tree_pt = (*it);
6200 deleted_element_pt.push_back(tree_pt->
object_pt());
6212 for (
int domain = 0; domain < n_proc; domain++)
6214 unsigned nelem = tmp_root_halo_element_pt[domain].size();
6215 for (
unsigned e = 0;
e < nelem;
e++)
6218 tmp_root_halo_element_pt[domain][
e]);
6221 nelem = tmp_root_haloed_element_pt[domain].size();
6222 for (
unsigned e = 0;
e < nelem;
e++)
6225 tmp_root_haloed_element_pt[domain][
e]);
6237 for (
unsigned e = 0;
e < nelem;
e++)
6242 unsigned nnod = el_pt->
nnode();
6243 for (
unsigned j = 0; j < nnod; j++)
6255 nnod = backed_up_nod_pt.size();
6256 for (
unsigned j = 0; j < nnod; j++)
6258 Node* nod_pt = backed_up_nod_pt[j];
6286 if (ref_mesh_pt != 0)
6297 <<
"Time for local rebuild of mesh from retained elements in "
6298 <<
" Mesh::prune_halo_elements_and_nodes(): " << t_end - t_start
6309 oomph_info <<
"Time for Mesh::classify_halo_and_haloed_nodes() "
6310 <<
" Mesh::prune_halo_elements_and_nodes(): "
6311 << t_end - t_start << std::endl;
6323 <<
" " << doc_info.
number() << std::endl;
6337 oomph_info <<
"Time for Mesh::reorder_nodes() "
6338 <<
" Mesh::prune_halo_elements_and_nodes(): "
6339 << t_end - t_start << std::endl;
6349 <<
"\n----------------------------------------------------\n";
6350 oomph_info <<
"After pruning: Processor " << my_rank <<
" holds "
6351 << this->
nelement() <<
" elements of which "
6353 <<
" are root halo elements \n while "
6355 <<
" are root haloed elements" << std::endl;
6358 oomph_info <<
"After pruning: Processor " << my_rank <<
" holds "
6362 <<
" are shared nodes." << std::endl;
6366 for (
int iproc = 0; iproc < n_proc; iproc++)
6368 oomph_info <<
"After pruning: With process " << iproc
6370 <<
" halo nodes, and " << std::endl
6376 <<
"----------------------------------------------------\n\n";
6391 double& max_efficiency,
6392 double& min_efficiency)
6395 int n_proc =
Comm_pt->nproc();
6396 int my_rank =
Comm_pt->my_rank();
6404 for (
int d = 0; d < n_proc; d++)
6407 count += halo_elem_pt.size();
6411 int nhalo_element_local = count;
6418 MPI_Gather(&nhalo_element_local,
6426 MPI_Gather(&n_element_local,
6436 av_efficiency = 0.0;
6438 double min = 1000000000.0;
6442 for (
int i = 0;
i < n_proc;
i++)
6444 unsigned nel = n_elements[
i];
6448 eff = double(n_elements[
i] - nhalo_elements[
i]) / double(nel);
6450 av_efficiency += eff;
6451 if (eff > max) max = eff;
6452 if (eff < min) min = eff;
6454 av_efficiency /= double(n_proc);
6458 MPI_Bcast(&max, 1, MPI_DOUBLE, 0,
Comm_pt->mpi_comm());
6459 MPI_Bcast(&min, 1, MPI_DOUBLE, 0,
Comm_pt->mpi_comm());
6460 MPI_Bcast(&av_efficiency, 1, MPI_DOUBLE, 0,
Comm_pt->mpi_comm());
6462 max_efficiency = max;
6463 min_efficiency = min;
6473 int my_rank =
Comm_pt->my_rank();
6474 int n_proc =
Comm_pt->nproc();
6476 std::ostringstream filename;
6477 std::ostringstream filename2;
6478 std::ofstream some_file;
6479 std::ofstream some_file2;
6483 <<
"elements_on_proc" << my_rank <<
"_" << doc_info.
number()
6485 some_file.open(filename.str().c_str());
6486 this->
output(some_file, 5);
6492 <<
"non_halo_elements_on_proc" << my_rank <<
"_"
6493 << doc_info.
number() <<
".dat";
6494 some_file.open(filename.str().c_str());
6498 for (
unsigned e = 0;
e < nelem;
e++)
6509 some_file <<
"Generalised Element " <<
e <<
"\n";
6514 f_el_pt->
output(some_file, 5);
6525 <<
"halo_elements_on_proc" << my_rank <<
"_" << doc_info.
number()
6527 some_file.open(filename.str().c_str());
6528 for (
int domain = 0; domain < n_proc; domain++)
6532 <<
"halo_elements_with_proc" << domain <<
"_on_proc" << my_rank
6533 <<
"_" << doc_info.
number() <<
".dat";
6534 some_file2.open(filename2.str().c_str());
6538 unsigned nelem = halo_elem_pt.size();
6539 for (
unsigned e = 0;
e < nelem;
e++)
6553 std::ostringstream error_message;
6555 <<
"Haloed element is not a leaf element. This shouldn't happen"
6557 error_message <<
"Here are the nodal positions: " << std::endl;
6558 unsigned nnod = ref_el_pt->
nnode();
6559 for (
unsigned j = 0; j < nnod; j++)
6562 unsigned n_dim = nod_pt->
ndim();
6563 for (
unsigned i = 0;
i < n_dim;
i++)
6565 error_message << nod_pt->
x(
i) <<
" ";
6567 error_message <<
"\n";
6569 OOMPH_CURRENT_FUNCTION,
6570 OOMPH_EXCEPTION_LOCATION);
6576 f_el_pt->
output(some_file, 5);
6577 f_el_pt->
output(some_file2, 5);
6582 some_file <<
"Generalised Element " <<
e <<
"\n";
6583 some_file2 <<
"Generalised Element " <<
e <<
"\n";
6593 <<
"haloed_elements_on_proc" << my_rank <<
"_" << doc_info.
number()
6595 some_file.open(filename.str().c_str());
6596 for (
int domain = 0; domain < n_proc; domain++)
6600 <<
"haloed_elements_with_proc" << domain <<
"_on_proc"
6601 << my_rank <<
"_" << doc_info.
number() <<
".dat";
6602 some_file2.open(filename2.str().c_str());
6607 unsigned nelem = haloed_elem_pt.size();
6608 for (
unsigned e = 0;
e < nelem;
e++)
6613 if (finite_el_pt != 0)
6624 std::ostringstream error_message;
6626 <<
"Haloed element is not a leaf element. This shouldn't happen"
6628 error_message <<
"Here are the nodal positions: " << std::endl;
6629 unsigned nnod = ref_el_pt->
nnode();
6630 for (
unsigned j = 0; j < nnod; j++)
6633 unsigned n_dim = nod_pt->
ndim();
6634 for (
unsigned i = 0;
i < n_dim;
i++)
6636 error_message << nod_pt->
x(
i) <<
" ";
6638 error_message <<
"\n";
6640 OOMPH_CURRENT_FUNCTION,
6641 OOMPH_EXCEPTION_LOCATION);
6647 finite_el_pt->
output(some_file, 5);
6648 finite_el_pt->
output(some_file2, 5);
6653 some_file <<
"Generalised Element " <<
e <<
"\n";
6654 some_file2 <<
"Generalised Element " <<
e <<
"\n";
6665 <<
"non_halo_nodes_on_proc" << my_rank <<
"_" << doc_info.
number()
6667 some_file.open(filename.str().c_str());
6668 unsigned nnod = this->
nnode();
6669 for (
unsigned j = 0; j < nnod; j++)
6674 unsigned ndim = nod_pt->
ndim();
6675 for (
unsigned i = 0;
i < ndim;
i++)
6677 some_file << nod_pt->
x(
i) <<
" ";
6690 <<
"nodes_on_proc" << my_rank <<
"_" << doc_info.
number()
6692 some_file.open(filename.str().c_str());
6693 for (
unsigned j = 0; j < nnod; j++)
6696 unsigned ndim = nod_pt->
ndim();
6697 for (
unsigned i = 0;
i < ndim;
i++)
6699 some_file << nod_pt->
x(
i) <<
" ";
6709 <<
"solid_nodes_on_proc" << my_rank <<
"_" << doc_info.
number()
6711 some_file.open(filename.str().c_str());
6712 unsigned nsnod = this->
nnode();
6713 for (
unsigned j = 0; j < nsnod; j++)
6717 if (solid_nod_pt != 0)
6719 unsigned ndim = solid_nod_pt->
ndim();
6720 for (
unsigned i = 0;
i < ndim;
i++)
6722 some_file << nod_pt->
x(
i) <<
" ";
6735 <<
"halo_nodes_on_proc" << my_rank <<
"_" << doc_info.
number()
6737 some_file.open(filename.str().c_str());
6738 for (
int domain = 0; domain < n_proc; domain++)
6742 <<
"halo_nodes_with_proc" << domain <<
"_on_proc" << my_rank
6743 <<
"_" << doc_info.
number() <<
".dat";
6744 some_file2.open(filename2.str().c_str());
6746 for (
unsigned j = 0; j < nnod; j++)
6749 unsigned ndim = nod_pt->
ndim();
6750 for (
unsigned i = 0;
i < ndim;
i++)
6752 some_file << nod_pt->
x(
i) <<
" ";
6753 some_file2 << nod_pt->
x(
i) <<
" ";
6770 <<
"haloed_nodes_on_proc" << my_rank <<
"_" << doc_info.
number()
6772 some_file.open(filename.str().c_str());
6773 for (
int domain = 0; domain < n_proc; domain++)
6777 <<
"haloed_nodes_with_proc" << domain <<
"_on_proc" << my_rank
6778 <<
"_" << doc_info.
number() <<
".dat";
6779 some_file2.open(filename2.str().c_str());
6781 for (
unsigned j = 0; j < nnod; j++)
6784 unsigned ndim = nod_pt->
ndim();
6785 for (
unsigned i = 0;
i < ndim;
i++)
6787 some_file << nod_pt->
x(
i) <<
" ";
6788 some_file2 << nod_pt->
x(
i) <<
" ";
6805 <<
"shared_nodes_on_proc" << my_rank <<
"_" << doc_info.
number()
6807 some_file.open(filename.str().c_str());
6808 for (
int domain = 0; domain < n_proc; domain++)
6812 <<
"shared_nodes_with_proc" << domain <<
"_on_proc" << my_rank
6813 <<
"_" << doc_info.
number() <<
".dat";
6814 some_file2.open(filename2.str().c_str());
6816 for (
unsigned j = 0; j < nnod; j++)
6819 unsigned ndim = nod_pt->
ndim();
6820 for (
unsigned i = 0;
i < ndim;
i++)
6822 some_file << nod_pt->
x(
i) <<
" ";
6823 some_file2 << nod_pt->
x(
i) <<
" ";
6826 <<
" " << nod_pt->
is_halo() <<
"\n";
6838 filename << doc_info.
directory() <<
"/" << doc_info.
label() <<
"mesh"
6839 << my_rank <<
"_" << doc_info.
number() <<
".dat";
6840 some_file.open(filename.str().c_str());
6841 this->
output(some_file, 5);
6847 filename << doc_info.
directory() <<
"/" << doc_info.
label() <<
"boundaries"
6848 << my_rank <<
"_" << doc_info.
number() <<
".dat";
6849 some_file.open(filename.str().c_str());
6860 for (
unsigned b = 0; b < nbound; b++)
6864 <<
"boundary_elements" << my_rank <<
"_" << b <<
"_"
6865 << doc_info.
number() <<
".dat";
6866 some_file.open(filename.str().c_str());
6868 for (
unsigned e = 0;
e < nelem;
e++)
6882 double& max_permitted_error_for_halo_check)
6884 oomph_info <<
"Doing check_halo_schemes for mesh: " <<
typeid(*this).name()
6892 std::ostringstream filename;
6893 std::ofstream shared_file;
6894 std::ofstream halo_file;
6895 std::ofstream haloed_file;
6896 std::ofstream ext_halo_file;
6897 std::ofstream ext_haloed_file;
6900 int n_proc =
Comm_pt->nproc();
6901 int my_rank =
Comm_pt->my_rank();
6912 for (
int dd = 0; dd < n_proc; dd++)
6915 filename << doc_info.
directory() <<
"/shared_node_check"
6916 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
6917 << dd <<
"_" << doc_info.
number() <<
".dat";
6918 shared_file.open(filename.str().c_str());
6919 shared_file <<
"ZONE " << std::endl;
6922 for (
unsigned j = 0; j < nnod; j++)
6925 unsigned ndim = nod_pt->
ndim();
6926 for (
unsigned i = 0;
i < ndim;
i++)
6928 shared_file << nod_pt->
position(
i) <<
" ";
6930 shared_file << j <<
" " << nod_pt << std::endl;
6934 if ((nnod == 0) && (
nelement() != 0))
6940 shared_file <<
"0.0" << std::endl;
6947 shared_file <<
" 1.0 1.1 " << std::endl;
6951 shared_file <<
" 1.0 1.1 1.1" << std::endl;
6955 shared_file.close();
6961 for (
int d = 0; d < n_proc; d++)
6964 std::set<Node*> shared_node_set;
6965 for (
unsigned i = 0;
i < n_vector;
i++)
6967 unsigned old_size = shared_node_set.size();
6969 unsigned new_size = shared_node_set.size();
6970 if (old_size == new_size)
6973 oomph_info <<
"Repeated node in shared node lookup scheme: " <<
i
6974 <<
"-th node with proc " << d <<
" : " << nod_pt
6976 unsigned n = nod_pt->
ndim();
6977 for (
unsigned ii = 0; ii < n; ii++)
6984 unsigned n_set = shared_node_set.size();
6985 if (n_vector != n_set)
6987 std::ostringstream warning_stream;
6989 <<
"WARNING: " << std::endl
6990 <<
"There seem to be repeated entries in shared node scheme "
6991 <<
"with proc " << d <<
".\n"
6992 <<
"n_vector=" << n_vector <<
"; n_set=" << n_set << std::endl;
6995 warning_stream <<
"Shared node scheme has been output in files like\n"
6996 << filename.str() << std::endl;
7000 warning_stream <<
"Re-run with doc_info enabled to see where the "
7001 "shared nodes are.\n";
7004 "Mesh::check_halo_schemes",
7005 OOMPH_EXCEPTION_LOCATION);
7012 double max_error = 0.0;
7015 for (
int d = 0; d < n_proc; d++)
7021 for (
int dd = 0; dd < n_proc; dd++)
7029 if (nnod_shared != 0)
7035 &nnod_share, 1, MPI_INT, dd, 0,
Comm_pt->mpi_comm(), &status);
7037 if (nnod_shared != nnod_share)
7039 std::ostringstream error_message;
7041 error_message <<
"Clash in numbers of shared nodes! "
7043 error_message <<
"# of shared nodes on proc " << dd <<
": "
7044 << nnod_shared << std::endl;
7045 error_message <<
"# of shared nodes on proc " << d <<
": "
7046 << nnod_share << std::endl;
7047 error_message <<
"(Re-)run Problem::check_halo_schemes() with "
7050 error_message <<
"to identify the problem" << std::endl;
7052 OOMPH_CURRENT_FUNCTION,
7053 OOMPH_EXCEPTION_LOCATION);
7061 MPI_Recv(&other_nodal_positions[0],
7062 nod_dim * nnod_share,
7071 for (
int j = 0; j < nnod_share; j++)
7074 for (
unsigned i = 0;
i < nod_dim;
i++)
7079 for (
unsigned i = 0;
i < nod_dim;
i++)
7081 x_share[
i] = other_nodal_positions[count];
7085 for (
unsigned i = 0;
i < nod_dim;
i++)
7088 (x_shared[
i] - x_share[
i]) * (x_shared[
i] - x_share[
i]);
7090 error = sqrt(error);
7093 if (error > max_permitted_error_for_halo_check)
7096 for (
unsigned i = 0;
i < nod_dim;
i++)
7100 for (
unsigned i = 0;
i < nod_dim;
i++)
7110 <<
" is hanging with masters" << std::endl;
7111 for (
unsigned m = 0;
7125 if (error > max_error)
7140 if (nnod_share != 0)
7144 MPI_Send(&nnod_share, 1, MPI_INT, d, 0,
Comm_pt->mpi_comm());
7151 for (
int j = 0; j < nnod_share; j++)
7153 for (
unsigned i = 0;
i < nod_dim;
i++)
7161 MPI_Send(&nodal_positions[0],
7162 nod_dim * nnod_share,
7171 oomph_info <<
"Max. error for shared nodes " << max_error << std::endl;
7172 if (max_error > max_permitted_error_for_halo_check)
7174 std::ostringstream error_message;
7175 error_message <<
"This is bigger than the permitted threshold "
7176 << max_permitted_error_for_halo_check << std::endl;
7178 <<
"If you believe this to be acceptable for your problem\n"
7179 <<
"increase Problem::Max_permitted_error_for_halo_check and re-run \n";
7181 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
7191 for (
int dd = 0; dd < n_proc; dd++)
7194 filename << doc_info.
directory() <<
"/halo_element_check"
7195 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
7196 << dd <<
"_" << doc_info.
number() <<
".dat";
7197 halo_file.open(filename.str().c_str());
7202 unsigned nelem = halo_elem_pt.size();
7204 for (
unsigned e = 0;
e < nelem;
e++)
7206 halo_file <<
"ZONE " << std::endl;
7210 if (finite_el_pt != 0)
7212 unsigned nnod = finite_el_pt->
nnode();
7213 for (
unsigned j = 0; j < nnod; j++)
7216 unsigned ndim = nod_pt->
ndim();
7217 for (
unsigned i = 0;
i < ndim;
i++)
7219 halo_file << nod_pt->
position(
i) <<
" ";
7221 halo_file << std::endl;
7229 for (
int d = 0; d < n_proc; d++)
7232 filename << doc_info.
directory() <<
"/haloed_element_check"
7233 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
7234 << d <<
"_" << doc_info.
number() <<
".dat";
7235 haloed_file.open(filename.str().c_str());
7240 unsigned nelem2 = haloed_elem_pt.size();
7241 for (
unsigned e = 0;
e < nelem2;
e++)
7243 haloed_file <<
"ZONE " << std::endl;
7247 if (finite_el_pt != 0)
7257 std::ostringstream error_message;
7258 error_message <<
"Haloed element is not a leaf element. This "
7261 error_message <<
"Here are the nodal positions: " << std::endl;
7262 unsigned nnod = ref_el_pt->
nnode();
7263 for (
unsigned j = 0; j < nnod; j++)
7266 unsigned n_dim = nod_pt->
ndim();
7267 for (
unsigned i = 0;
i < n_dim;
i++)
7269 error_message << nod_pt->
x(
i) <<
" ";
7271 error_message <<
"\n";
7273 OOMPH_CURRENT_FUNCTION,
7274 OOMPH_EXCEPTION_LOCATION);
7278 unsigned nnod2 = finite_el_pt->
nnode();
7279 for (
unsigned j = 0; j < nnod2; j++)
7282 unsigned ndim = nod_pt->
ndim();
7283 for (
unsigned i = 0;
i < ndim;
i++)
7285 haloed_file << nod_pt->
position(
i) <<
" ";
7287 haloed_file << std::endl;
7291 haloed_file.close();
7299 bool shout_and_terminate =
false;
7302 for (
int d = 0; d < n_proc; d++)
7308 for (
int dd = 0; dd < n_proc; dd++)
7318 int nelem_haloed = haloed_elem_pt.size();
7320 if (nelem_haloed != 0)
7326 &nelem_halo, 1, MPI_INT, dd, 0,
Comm_pt->mpi_comm(), &status);
7327 if (nelem_halo != nelem_haloed)
7329 std::ostringstream error_message;
7331 <<
"Clash in numbers of halo and haloed elements! "
7333 error_message <<
"# of haloed elements whose halo counterpart "
7335 << dd <<
": " << nelem_haloed << std::endl;
7336 error_message <<
"# of halo elements whose non-halo "
7337 "counterpart lives on proc "
7338 << d <<
": " << nelem_halo << std::endl;
7339 error_message <<
"(Re-)run Problem::check_halo_schemes() with "
7342 error_message <<
"to identify the problem" << std::endl;
7344 OOMPH_CURRENT_FUNCTION,
7345 OOMPH_EXCEPTION_LOCATION);
7360 unsigned n_nodal_positions = 0;
7361 MPI_Recv(&n_nodal_positions,
7369 if (n_nodal_positions > 0)
7371 MPI_Recv(&other_nodal_positions[0],
7383 unsigned n_other_nodal_hangings = 0;
7384 MPI_Recv(&n_other_nodal_hangings,
7391 if (n_other_nodal_hangings > 0)
7393 other_nodal_hangings.resize(n_other_nodal_hangings);
7394 MPI_Recv(&other_nodal_hangings[0],
7395 n_other_nodal_hangings,
7407 filename << doc_info.
directory() <<
"/error_haloed_check"
7408 << doc_info.
label() <<
"_on_proc" << my_rank
7409 <<
"_with_proc" << dd <<
"_" << doc_info.
number()
7411 haloed_file.open(filename.str().c_str());
7413 filename << doc_info.
directory() <<
"/error_halo_check"
7414 << doc_info.
label() <<
"_on_proc" << my_rank
7415 <<
"_with_proc" << dd <<
"_" << doc_info.
number()
7417 halo_file.open(filename.str().c_str());
7421 unsigned count_hanging = 0;
7422 for (
int e = 0;
e < nelem_haloed;
e++)
7427 if (finite_el_pt != 0)
7429 unsigned nnod_this_el = finite_el_pt->
nnode();
7430 for (
unsigned j = 0; j < nnod_this_el; j++)
7438 for (
unsigned i = 0;
i < nod_dim;
i++)
7444 for (
unsigned i = 0;
i < nod_dim;
i++)
7446 x_halo[
i] = other_nodal_positions[count];
7452 for (
unsigned i = 0;
i < nod_dim;
i++)
7455 (x_haloed[
i] - x_halo[
i]) * (x_haloed[
i] - x_halo[
i]);
7457 error = sqrt(error);
7459 if (error > max_error)
7463 double tol = 1.0e-12;
7467 <<
"Discrepancy between nodal coordinates of halo(ed)"
7468 <<
"element larger than tolerance (" << tol
7469 <<
")\n Error: " << error <<
"\n";
7473 unsigned nval = nod_pt->
nvalue();
7474 int nval_other = other_nodal_hangings[count_hanging];
7476 if (
int(nval) != nval_other)
7479 <<
"Number of values of node, " << nval
7480 <<
", does not match number of values on other proc, "
7481 << nval_other << std::endl;
7483 shout_and_terminate =
true;
7487 int other_geom_hanging = 0;
7490 for (
int i = -1;
i < int(nval);
i++)
7492 int nmaster_other = other_nodal_hangings[count_hanging];
7496 if (
i == -1) other_geom_hanging = nmaster_other;
7503 if (
int(nmaster) != nmaster_other)
7506 <<
"Number of master nodes for hanging value "
7507 <<
i <<
" of node, " << nmaster
7508 <<
", does not match number of master "
7509 <<
"nodes on other proc, " << nmaster_other
7512 shout_and_terminate =
true;
7520 if (nmaster_other != 0)
7524 <<
" of node is not hanging whereas "
7525 <<
" node on other proc has " << nmaster_other
7526 <<
" masters and therefore is hanging. \n";
7528 shout_and_terminate =
true;
7539 <<
"Error(s) displayed above are for "
7540 <<
"domain with non-halo (i.e. haloed) elem: " << dd
7543 <<
"Domain with halo elem: " << d
7549 <<
"Current processor is " << my_rank <<
"\n"
7550 <<
"Nodal positions: " << x_halo[0] <<
"\n"
7551 <<
"and haloed: " << x_haloed[0]
7558 <<
"Current processor is " << my_rank <<
"\n"
7559 <<
"Nodal positions: " << x_halo[0] <<
" "
7560 << x_halo[1] <<
"\n"
7561 <<
"and haloed: " << x_haloed[0] <<
" "
7569 <<
"Current processor is " << my_rank <<
"\n"
7570 <<
"Nodal positions: " << x_halo[0] <<
" "
7571 << x_halo[1] <<
" " << x_halo[2] <<
"\n"
7572 <<
"and haloed: " << x_haloed[0] <<
" "
7573 << x_haloed[1] <<
" " << x_haloed[2]
7580 "Nodal dimension not equal to 1, 2 or 3\n",
7581 OOMPH_CURRENT_FUNCTION,
7582 OOMPH_EXCEPTION_LOCATION);
7589 for (
unsigned i = 0;
i < nod_dim;
i++)
7591 haloed_file << x_haloed[
i] <<
" ";
7592 halo_file << x_halo[
i] <<
" ";
7594 haloed_file << error <<
" " << my_rank <<
" " << dd
7598 halo_file << error <<
" " << my_rank <<
" " << dd
7599 <<
" " << other_geom_hanging << std::endl;
7609 haloed_file.close();
7626 unsigned nelem_halo = halo_elem_pt.size();
7628 if (nelem_halo != 0)
7632 MPI_Send(&nelem_halo, 1, MPI_UNSIGNED, d, 0,
Comm_pt->mpi_comm());
7645 nodal_positions.reserve(nod_dim * nnod_first_el * nelem_halo);
7651 for (
unsigned e = 0;
e < nelem_halo;
e++)
7655 if (finite_el_pt != 0)
7657 unsigned nnod_this_el = finite_el_pt->
nnode();
7658 for (
unsigned j = 0; j < nnod_this_el; j++)
7667 oomph_info <<
"element: " << finite_el_pt << std::endl;
7668 for (
unsigned i = 0;
i < finite_el_pt->
nnode();
i++)
7675 OOMPH_CURRENT_FUNCTION,
7676 OOMPH_EXCEPTION_LOCATION);
7681 for (
unsigned i = 0;
i < nod_dim;
i++)
7685 nodal_positions.push_back(nod_pt->
position(
i));
7688 unsigned nval = nod_pt->
nvalue();
7689 nodal_hangings.push_back(nval);
7690 for (
int i = -1;
i < int(nval);
i++)
7695 nodal_hangings.push_back(nmaster);
7699 nodal_hangings.push_back(0);
7707 unsigned n_nodal_positions = nodal_positions.size();
7710 unsigned n_nodal_hangings = nodal_hangings.size();
7717 &n_nodal_positions, 1, MPI_UNSIGNED, d, 0,
Comm_pt->mpi_comm());
7718 if (n_nodal_positions > 0)
7720 MPI_Send(&nodal_positions[0],
7728 &n_nodal_hangings, 1, MPI_UNSIGNED, d, 1,
Comm_pt->mpi_comm());
7729 if (n_nodal_hangings > 0)
7731 MPI_Send(&nodal_hangings[0],
7743 oomph_info <<
"Max. error for halo/haloed elements " << max_error
7745 if (max_error > max_permitted_error_for_halo_check)
7747 shout_and_terminate =
true;
7748 oomph_info <<
"This is bigger than the permitted threshold "
7749 << max_permitted_error_for_halo_check << std::endl;
7751 <<
"If you believe this to be acceptable for your problem\n"
7752 <<
"increase Problem::Max_permitted_error_for_halo_check and re-run \n";
7755 if (shout_and_terminate)
7758 OOMPH_CURRENT_FUNCTION,
7759 OOMPH_EXCEPTION_LOCATION);
7769 for (
int dd = 0; dd < n_proc; dd++)
7772 filename << doc_info.
directory() <<
"/halo_node_check"
7773 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
7774 << dd <<
"_" << doc_info.
number() <<
".dat";
7775 halo_file.open(filename.str().c_str());
7776 halo_file <<
"ZONE " << std::endl;
7779 for (
unsigned j = 0; j < nnod; j++)
7782 unsigned ndim = nod_pt->
ndim();
7783 for (
unsigned i = 0;
i < ndim;
i++)
7785 halo_file << nod_pt->
position(
i) <<
" ";
7787 halo_file << nod_pt->
is_hanging() << std::endl;
7792 if ((nnod == 0) && (
nelement() != 0))
7798 halo_file <<
"0.0" << std::endl;
7805 halo_file <<
" 1.0 1.1 " << std::endl;
7809 halo_file <<
" 1.0 1.1 1.1" << std::endl;
7818 for (
int d = 0; d < n_proc; d++)
7821 filename << doc_info.
directory() <<
"/haloed_node_check"
7822 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
7823 << d <<
"_" << doc_info.
number() <<
".dat";
7824 haloed_file.open(filename.str().c_str());
7825 haloed_file <<
"ZONE " << std::endl;
7828 for (
unsigned j = 0; j < nnod; j++)
7831 unsigned ndim = nod_pt->
ndim();
7832 for (
unsigned i = 0;
i < ndim;
i++)
7834 haloed_file << nod_pt->
position(
i) <<
" ";
7836 haloed_file << nod_pt->
is_hanging() << std::endl;
7840 if ((nnod == 0) && (
nelement() != 0))
7846 haloed_file <<
"0.0" << std::endl;
7853 haloed_file <<
" 1.0 1.1 " << std::endl;
7857 haloed_file <<
" 1.0 1.1 1.1" << std::endl;
7861 haloed_file.close();
7870 for (
int d = 0; d < n_proc; d++)
7876 for (
int dd = 0; dd < n_proc; dd++)
7885 if (nnod_haloed != 0)
7891 &nnod_halo, 1, MPI_INT, dd, 0,
Comm_pt->mpi_comm(), &status);
7893 if (nnod_haloed != nnod_halo)
7895 std::ostringstream error_message;
7897 error_message <<
"Clash in numbers of halo and haloed nodes! "
7900 <<
"# of haloed nodes whose halo counterpart lives on proc "
7901 << dd <<
": " << nnod_haloed << std::endl;
7903 <<
"# of halo nodes whose non-halo counterpart lives on proc "
7904 << d <<
": " << nnod_halo << std::endl;
7906 <<
"(Re-)run Mesh::check_halo_schemes() with DocInfo object"
7908 error_message <<
"to identify the problem" << std::endl;
7910 OOMPH_CURRENT_FUNCTION,
7911 OOMPH_EXCEPTION_LOCATION);
7919 MPI_Recv(&other_nodal_positions[0],
7920 nod_dim * nnod_halo,
7929 for (
int j = 0; j < nnod_halo; j++)
7932 for (
unsigned i = 0;
i < nod_dim;
i++)
7937 for (
unsigned i = 0;
i < nod_dim;
i++)
7939 x_halo[
i] = other_nodal_positions[count];
7943 for (
unsigned i = 0;
i < nod_dim;
i++)
7946 (x_haloed[
i] - x_halo[
i]) * (x_haloed[
i] - x_halo[
i]);
7948 error = sqrt(error);
7949 if (error > max_error)
7968 MPI_Send(&nnod_halo, 1, MPI_INT, d, 0,
Comm_pt->mpi_comm());
7975 for (
int j = 0; j < nnod_halo; j++)
7977 for (
unsigned i = 0;
i < nod_dim;
i++)
7985 MPI_Send(&nodal_positions[0],
7986 nod_dim * nnod_halo,
7995 oomph_info <<
"Max. error for halo/haloed nodes " << max_error << std::endl;
7997 if (max_error > max_permitted_error_for_halo_check)
7999 std::ostringstream error_message;
8000 error_message <<
"This is bigger than the permitted threshold "
8001 << max_permitted_error_for_halo_check << std::endl;
8003 <<
"If you believe this to be acceptable for your problem\n"
8004 <<
"increase Problem::Max_permitted_error_for_halo_check and re-run \n";
8006 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
8017 for (
int dd = 0; dd < n_proc; dd++)
8020 filename << doc_info.
directory() <<
"/ext_halo_element_check"
8021 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
8022 << dd <<
"_" << doc_info.
number() <<
".dat";
8023 ext_halo_file.open(filename.str().c_str());
8025 ext_halo_file.close();
8029 filename << doc_info.
directory() <<
"/ext_halo_node_check"
8030 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
8031 << dd <<
"_" << doc_info.
number() <<
".dat";
8032 ext_halo_file.open(filename.str().c_str());
8038 unsigned nelem = ext_halo_elem_pt.size();
8040 for (
unsigned e = 0;
e < nelem;
e++)
8042 ext_halo_file <<
"ZONE " << std::endl;
8046 if (finite_el_pt != 0)
8048 unsigned nnod = finite_el_pt->
nnode();
8049 for (
unsigned j = 0; j < nnod; j++)
8052 unsigned ndim = nod_pt->
ndim();
8053 for (
unsigned i = 0;
i < ndim;
i++)
8055 ext_halo_file << nod_pt->
position(
i) <<
" ";
8057 ext_halo_file << std::endl;
8061 ext_halo_file.close();
8065 for (
int d = 0; d < n_proc; d++)
8068 filename << doc_info.
directory() <<
"/ext_haloed_element_check"
8069 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
8070 << d <<
"_" << doc_info.
number() <<
".dat";
8071 ext_haloed_file.open(filename.str().c_str());
8073 ext_haloed_file.close();
8077 filename << doc_info.
directory() <<
"/ext_haloed_node_check"
8078 << doc_info.
label() <<
"_on_proc" << my_rank <<
"_with_proc"
8079 << d <<
"_" << doc_info.
number() <<
".dat";
8080 ext_haloed_file.open(filename.str().c_str());
8086 unsigned nelem2 = ext_haloed_elem_pt.size();
8087 for (
unsigned e = 0;
e < nelem2;
e++)
8089 ext_haloed_file <<
"ZONE " << std::endl;
8093 if (finite_el_pt != 0)
8095 unsigned nnod2 = finite_el_pt->
nnode();
8096 for (
unsigned j = 0; j < nnod2; j++)
8099 unsigned ndim = nod_pt->
ndim();
8100 for (
unsigned i = 0;
i < ndim;
i++)
8102 ext_haloed_file << nod_pt->
position(
i) <<
" ";
8104 ext_haloed_file << std::endl;
8108 ext_haloed_file.close();
8116 shout_and_terminate =
false;
8119 for (
int d = 0; d < n_proc; d++)
8125 for (
int dd = 0; dd < n_proc; dd++)
8136 int nelem_haloed = ext_haloed_elem_pt.size();
8138 if (nelem_haloed != 0)
8144 &nelem_halo, 1, MPI_INT, dd, 0,
Comm_pt->mpi_comm(), &status);
8145 if (nelem_halo != nelem_haloed)
8147 std::ostringstream error_message;
8149 <<
"Clash in numbers of external halo and haloed elements! "
8151 error_message <<
"# of external haloed elements whose halo "
8152 "counterpart lives on proc "
8153 << dd <<
": " << nelem_haloed << std::endl;
8154 error_message <<
"# of external halo elements whose non-halo "
8155 "counterpart lives on proc "
8156 << d <<
": " << nelem_halo << std::endl;
8157 error_message <<
"(Re-)run Problem::check_halo_schemes() with "
8160 error_message <<
"to identify the problem" << std::endl;
8162 OOMPH_CURRENT_FUNCTION,
8163 OOMPH_EXCEPTION_LOCATION);
8180 unsigned n_nodal_positions = 0;
8181 MPI_Recv(&n_nodal_positions,
8189 if (n_nodal_positions > 0)
8191 MPI_Recv(&other_nodal_positions[0],
8203 unsigned n_other_nodal_hangings = 0;
8204 MPI_Recv(&n_other_nodal_hangings,
8211 if (n_other_nodal_hangings > 0)
8213 other_nodal_hangings.resize(n_other_nodal_hangings);
8214 MPI_Recv(&other_nodal_hangings[0],
8215 n_other_nodal_hangings,
8227 filename << doc_info.
directory() <<
"/error_ext_haloed_check"
8228 << doc_info.
label() <<
"_on_proc" << my_rank
8229 <<
"_with_proc" << dd <<
"_" << doc_info.
number()
8231 ext_haloed_file.open(filename.str().c_str());
8233 filename << doc_info.
directory() <<
"/error_ext_halo_check"
8234 << doc_info.
label() <<
"_on_proc" << my_rank
8235 <<
"_with_proc" << dd <<
"_" << doc_info.
number()
8237 ext_halo_file.open(filename.str().c_str());
8241 unsigned count_hanging = 0;
8242 for (
int e = 0;
e < nelem_haloed;
e++)
8247 if (finite_el_pt != 0)
8249 unsigned nnod_this_el = finite_el_pt->
nnode();
8250 for (
unsigned j = 0; j < nnod_this_el; j++)
8258 for (
unsigned i = 0;
i < nod_dim;
i++)
8264 for (
unsigned i = 0;
i < nod_dim;
i++)
8266 x_halo[
i] = other_nodal_positions[count];
8272 for (
unsigned i = 0;
i < nod_dim;
i++)
8275 (x_haloed[
i] - x_halo[
i]) * (x_haloed[
i] - x_halo[
i]);
8277 error = sqrt(error);
8279 if (error > max_error)
8283 double tol = 1.0e-12;
8286 oomph_info <<
"Discrepancy between nodal coordinates "
8287 "of external halo(ed)"
8288 <<
"element larger than tolerance (" << tol
8289 <<
")\n Error: " << error <<
"\n";
8293 unsigned nval = nod_pt->
nvalue();
8294 int nval_other = other_nodal_hangings[count_hanging];
8296 if (
int(nval) != nval_other)
8299 <<
"Number of values of node, " << nval
8300 <<
", does not match number of values on other proc, "
8301 << nval_other << std::endl;
8303 shout_and_terminate =
true;
8307 int other_geom_hanging = 0;
8310 for (
int i = -1;
i < int(nval);
i++)
8312 int nmaster_other = other_nodal_hangings[count_hanging];
8316 if (
i == -1) other_geom_hanging = nmaster_other;
8323 if (
int(nmaster) != nmaster_other)
8326 <<
"Number of master nodes for hanging value "
8327 <<
i <<
" of node, " << nmaster
8328 <<
", does not match number of master "
8329 <<
"nodes on other proc, " << nmaster_other
8332 shout_and_terminate =
true;
8340 if (nmaster_other != 0)
8344 <<
" of node is not hanging whereas "
8345 <<
" node on other proc has " << nmaster_other
8346 <<
" masters and therefore is hanging. \n";
8348 shout_and_terminate =
true;
8358 oomph_info <<
"Error(s) displayed above are for "
8359 <<
"domain with external non-halo (i.e. "
8363 <<
"Domain with halo elem: " << d
8369 <<
"Current processor is " << my_rank <<
"\n"
8370 <<
"Nodal positions: " << x_halo[0] <<
"\n"
8371 <<
"and haloed: " << x_haloed[0]
8378 <<
"Current processor is " << my_rank <<
"\n"
8379 <<
"Nodal positions: " << x_halo[0] <<
" "
8380 << x_halo[1] <<
"\n"
8381 <<
"and haloed: " << x_haloed[0] <<
" "
8389 <<
"Current processor is " << my_rank <<
"\n"
8390 <<
"Nodal positions: " << x_halo[0] <<
" "
8391 << x_halo[1] <<
" " << x_halo[2] <<
"\n"
8392 <<
"and haloed: " << x_haloed[0] <<
" "
8393 << x_haloed[1] <<
" " << x_haloed[2]
8400 "Nodal dimension not equal to 1, 2 or 3\n",
8401 OOMPH_CURRENT_FUNCTION,
8402 OOMPH_EXCEPTION_LOCATION);
8409 for (
unsigned i = 0;
i < nod_dim;
i++)
8411 ext_haloed_file << x_haloed[
i] <<
" ";
8412 ext_halo_file << x_halo[
i] <<
" ";
8415 << error <<
" " << my_rank <<
" " << dd <<
" "
8418 ext_halo_file << error <<
" " << my_rank <<
" " << dd
8419 <<
" " << other_geom_hanging
8430 ext_haloed_file.close();
8431 ext_halo_file.close();
8448 unsigned nelem_halo = ext_halo_elem_pt.size();
8450 if (nelem_halo != 0)
8454 MPI_Send(&nelem_halo, 1, MPI_UNSIGNED, d, 0,
Comm_pt->mpi_comm());
8462 unsigned nnod_first_el = fe_pt->
nnode();
8465 nodal_positions.reserve(nod_dim * nnod_first_el * nelem_halo);
8471 for (
unsigned e = 0;
e < nelem_halo;
e++)
8475 if (finite_el_pt != 0)
8477 unsigned nnod_this_el = finite_el_pt->
nnode();
8478 for (
unsigned j = 0; j < nnod_this_el; j++)
8484 for (
unsigned i = 0;
i < nod_dim;
i++)
8486 nodal_positions.push_back(nod_pt->
position(
i));
8490 unsigned nval = nod_pt->
nvalue();
8491 nodal_hangings.push_back(nval);
8492 for (
int i = -1;
i < int(nval);
i++)
8497 nodal_hangings.push_back(nmaster);
8501 nodal_hangings.push_back(0);
8509 unsigned n_nodal_positions = nodal_positions.size();
8512 unsigned n_nodal_hangings = nodal_hangings.size();
8519 &n_nodal_positions, 1, MPI_UNSIGNED, d, 0,
Comm_pt->mpi_comm());
8520 if (n_nodal_positions > 0)
8522 MPI_Send(&nodal_positions[0],
8530 &n_nodal_hangings, 1, MPI_UNSIGNED, d, 1,
Comm_pt->mpi_comm());
8531 if (n_nodal_hangings > 0)
8533 MPI_Send(&nodal_hangings[0],
8545 oomph_info <<
"Max. error for external halo/haloed elements " << max_error
8547 if (max_error > max_permitted_error_for_halo_check)
8549 shout_and_terminate =
true;
8550 oomph_info <<
"This is bigger than the permitted threshold "
8551 << max_permitted_error_for_halo_check << std::endl;
8553 <<
"If you believe this to be acceptable for your problem\n"
8554 <<
"increase Problem::Max_permitted_error_for_halo_check and re-run \n";
8557 if (shout_and_terminate)
8560 OOMPH_CURRENT_FUNCTION,
8561 OOMPH_EXCEPTION_LOCATION);
8573 unsigned n = ext_halo_node_pt.size();
8574 for (
unsigned j = 0; j < n; j++)
8576 if (ext_halo_node_pt[j] == nod_pt)
8592 int n_proc =
Comm_pt->nproc();
8593 int my_rank =
Comm_pt->my_rank();
8608 for (
int domain = 0; domain < n_proc; domain++)
8611 send_displacement[domain] = send_data.size();
8615 if (domain != my_rank)
8624 unsigned nnod = backup_pt.size();
8628 for (
unsigned j = 0; j < nnod; j++)
8631 Node* nod_pt = backup_pt[j];
8637 send_data.push_back(j);
8648 send_data.push_back(-1);
8651 send_n[domain] = send_data.size() - send_displacement[domain];
8660 &send_n[0], 1, MPI_INT, &receive_n[0], 1, MPI_INT,
Comm_pt->mpi_comm());
8666 int receive_data_count = 0;
8667 for (
int rank = 0; rank < n_proc; ++rank)
8670 receive_displacement[rank] = receive_data_count;
8671 receive_data_count += receive_n[rank];
8676 if (receive_data_count == 0)
8678 ++receive_data_count;
8684 if (send_data.size() == 0)
8686 send_data.resize(1);
8690 MPI_Alltoallv(&send_data[0],
8692 &send_displacement[0],
8696 &receive_displacement[0],
8701 for (
int send_rank = 0; send_rank < n_proc; send_rank++)
8705 if ((send_rank != my_rank) && (receive_n[send_rank] != 0))
8708 unsigned count = receive_displacement[send_rank];
8714 int next_one = receive_data[count++];
8734 unsigned nnod = backup_pt.size();
8738 for (
unsigned j = 0; j < nnod; j++)
8741 Node* nod_pt = backup_pt[j];
8767 int int_ndof_types = -1;
8771 int_ndof_types =
element_pt(0)->ndof_types();
8775 for (
unsigned i = 1;
i < nel;
i++)
8779 std::ostringstream error_message;
8781 <<
"Every element in the mesh must have the same number of "
8782 <<
"types of DOF for ndof_types() to work\n"
8783 <<
"Element 0 has " << int_ndof_types <<
" DOF types\n"
8784 <<
"Element " <<
i <<
" [out of a total of " << nel <<
" ] has "
8786 <<
"Element types are: Element 0:" <<
typeid(*
element_pt(0)).name()
8788 <<
" Current Element :" <<
typeid(*
element_pt(
i)).name()
8791 OOMPH_CURRENT_FUNCTION,
8792 OOMPH_EXCEPTION_LOCATION);
8798 #ifdef OOMPH_HAS_MPI
8809 unsigned nproc =
Comm_pt->nproc();
8810 unsigned my_rank =
Comm_pt->my_rank();
8815 int* ndof_types_recv = 0;
8818 ndof_types_recv =
new int[nproc];
8821 MPI_Gather(&int_ndof_types,
8837 for (
unsigned p = 1; p < nproc; p++)
8839 if (ndof_types_recv[p] != -1)
8844 if (int_ndof_types == -1)
8846 int_ndof_types = ndof_types_recv[p];
8850 else if (int_ndof_types != ndof_types_recv[p])
8852 std::ostringstream error_message;
8854 <<
"The elements in this mesh must have the same number "
8855 <<
"of types of DOF on each processor";
8856 for (
unsigned p = 0; p < nproc; p++)
8858 if (ndof_types_recv[p] != -1)
8860 error_message <<
"Processor " << p <<
" : "
8861 << ndof_types_recv[p] <<
"\n";
8865 error_message <<
"Processor " << p <<
" : (no elements)\n";
8869 OOMPH_CURRENT_FUNCTION,
8870 OOMPH_EXCEPTION_LOCATION);
8877 for (
unsigned p = 1; p < nproc; p++)
8879 if (ndof_types_recv[p] == -1)
8881 MPI_Send(&int_ndof_types, 1, MPI_INT, p, 0,
Comm_pt->mpi_comm());
8885 delete[] ndof_types_recv;
8890 else if (int_ndof_types == -1)
8892 MPI_Recv(&int_ndof_types,
8906 if (int_ndof_types == -1) int_ndof_types = 0;
8908 return unsigned(int_ndof_types);
8931 std::ostringstream error_message;
8933 <<
"Every element in the mesh must have the same number of "
8934 <<
"elemental dimension for elemental_dimension() to work.\n"
8935 <<
"Element 0 has elemental dimension " << int_dim <<
"\n"
8936 <<
"Element " <<
i <<
" has elemental dimension "
8939 OOMPH_CURRENT_FUNCTION,
8940 OOMPH_EXCEPTION_LOCATION);
8946 #ifdef OOMPH_HAS_MPI
8957 unsigned nproc =
Comm_pt->nproc();
8958 unsigned my_rank =
Comm_pt->my_rank();
8966 dim_recv =
new int[nproc];
8970 &int_dim, 1, MPI_INT, dim_recv, 1, MPI_INT, 0,
Comm_pt->mpi_comm());
8979 for (
unsigned p = 1; p < nproc; p++)
8981 if (dim_recv[p] != -1)
8988 int_dim = dim_recv[p];
8992 else if (int_dim != dim_recv[p])
8994 std::ostringstream error_message;
8996 <<
"The elements in this mesh must have the same elemental "
8997 <<
"dimension number on each processor";
8998 for (
unsigned p = 0; p < nproc; p++)
9000 if (dim_recv[p] != -1)
9002 error_message <<
"Processor " << p <<
" : " << dim_recv[p]
9007 error_message <<
"Processor " << p <<
" : (no elements)\n";
9011 OOMPH_CURRENT_FUNCTION,
9012 OOMPH_EXCEPTION_LOCATION);
9020 for (
unsigned p = 1; p < nproc; p++)
9022 if (dim_recv[p] == -1)
9024 MPI_Send(&int_dim, 1, MPI_INT, p, 0,
Comm_pt->mpi_comm());
9033 else if (int_dim == -1)
9036 &int_dim, 1, MPI_INT, 0, 0,
Comm_pt->mpi_comm(), MPI_STATUS_IGNORE);
9044 if (int_dim == -1) int_dim = 0;
9046 return unsigned(int_dim);
9069 std::ostringstream error_message;
9071 <<
"Every element in the mesh must have the same number of "
9072 <<
"nodal dimension for nodal_dimension() to work.\n"
9073 <<
"Element 0 has nodal dimension " << int_dim <<
"\n"
9074 <<
"Element " <<
i <<
" has nodal dimension "
9077 OOMPH_CURRENT_FUNCTION,
9078 OOMPH_EXCEPTION_LOCATION);
9084 #ifdef OOMPH_HAS_MPI
9095 unsigned nproc =
Comm_pt->nproc();
9096 unsigned my_rank =
Comm_pt->my_rank();
9104 dim_recv =
new int[nproc];
9108 &int_dim, 1, MPI_INT, dim_recv, 1, MPI_INT, 0,
Comm_pt->mpi_comm());
9117 for (
unsigned p = 1; p < nproc; p++)
9119 if (dim_recv[p] != -1)
9126 int_dim = dim_recv[p];
9130 else if (int_dim != dim_recv[p])
9132 std::ostringstream error_message;
9134 <<
"The elements in this mesh must have the same nodal "
9135 <<
"dimension number on each processor";
9136 for (
unsigned p = 0; p < nproc; p++)
9138 if (dim_recv[p] != -1)
9140 error_message <<
"Processor " << p <<
" : " << dim_recv[p]
9145 error_message <<
"Processor " << p <<
" : (no elements)\n";
9149 OOMPH_CURRENT_FUNCTION,
9150 OOMPH_EXCEPTION_LOCATION);
9158 for (
unsigned p = 1; p < nproc; p++)
9160 if (dim_recv[p] == -1)
9162 MPI_Send(&int_dim, 1, MPI_INT, p, 0,
Comm_pt->mpi_comm());
9171 else if (int_dim == -1)
9174 &int_dim, 1, MPI_INT, 0, 0,
Comm_pt->mpi_comm(), MPI_STATUS_IGNORE);
9182 if (int_dim == -1) int_dim = 0;
9184 return unsigned(int_dim);
9192 #ifdef OOMPH_HAS_MPI
9202 for (
unsigned i = 0;
i <
nnode();
i++)
9219 bool found_a_master_in_external_halo_storage =
false;
9220 for (
unsigned m = 0; m < hang_pt->
nmaster(); m++)
9226 bool found_this_master_in_external_halo_storage =
false;
9227 for (
int d = 0; d <
Comm_pt->nproc(); d++)
9238 found_this_master_in_external_halo_storage =
true;
9245 if (found_this_master_in_external_halo_storage)
9248 found_a_master_in_external_halo_storage =
true;
9255 if (found_a_master_in_external_halo_storage)
9265 unsigned n_value = nod_pt->
nvalue();
9267 unsigned n_dim = nod_pt->
ndim();
9270 for (
unsigned t = 0;
t < nt;
t++)
9272 nod_pt->
value(
t, values);
9273 for (
unsigned i = 0;
i < n_value;
i++)
9278 for (
unsigned i = 0;
i < n_dim;
i++)
9280 nod_pt->
x(
t,
i) = position[
i];
9287 if (alg_node_pt != 0)
9289 bool update_all_time_levels =
true;
9297 if (solid_node_pt != 0)
9299 unsigned n_lagrangian = solid_node_pt->
nlagrangian();
9300 for (
unsigned i = 0;
i < n_lagrangian;
i++)
9354 int d = (*it).first;
9358 for (
unsigned j = 0; j < n_ext_halo_nod; j++)
9362 for (
unsigned i_bnd = 0; i_bnd < n_bnd; i_bnd++)
9378 int d = (*it).first;
9381 for (
unsigned j = 0; j < n_ext_halo_nod; j++)
9384 bool is_a_mesh_node =
false;
9385 unsigned n_node =
nnode();
9386 for (
unsigned jj = 0; jj < n_node; jj++)
9390 is_a_mesh_node =
true;
9396 if (!is_a_mesh_node)
9407 int dd = (*itt).first;
9412 for (
unsigned jjj = 0; jjj < n_ext_halo; jjj++)
9417 is_a_mesh_node =
true;
9425 if (!is_a_mesh_node)
9439 int d = (*it).first;
9442 for (
unsigned e = 0;
e < n_ext_halo_el;
e++)
9460 #ifdef OOMPH_HAS_MPI
9482 bool already_external_haloed_element =
false;
9483 unsigned external_haloed_el_index = 0;
9484 for (
unsigned eh = 0; eh < n_extern_haloed; eh++)
9489 already_external_haloed_element =
true;
9491 external_haloed_el_index = eh;
9497 if (!already_external_haloed_element)
9502 return n_extern_haloed;
9507 return external_haloed_el_index;
9522 bool is_an_external_haloed_node =
false;
9523 unsigned external_haloed_node_index = 0;
9524 for (
unsigned k = 0; k < n_ext_haloed_nod; k++)
9528 is_an_external_haloed_node =
true;
9529 external_haloed_node_index = k;
9535 if (!is_an_external_haloed_node)
9540 return n_ext_haloed_nod;
9545 return external_haloed_node_index;
9567 unsigned long n_node =
nnode();
9570 for (
unsigned n = 0; n < n_node; n++)
9585 for (
unsigned k = 0; k < n_lagrangian_type; k++)
9588 for (
unsigned j = 0; j < n_lagrangian; j++)
9612 namespace ParaviewHelper
9617 pvd_file <<
"<?xml version=\"1.0\"?>" << std::endl
9618 <<
"<VTKFile type=\"Collection\" version=\"0.1\">" << std::endl
9619 <<
"<Collection>" << std::endl;
9629 pvd_file <<
"<DataSet timestep=\"" << time <<
"\" ";
9632 pvd_file <<
"part=\"0\" ";
9635 pvd_file <<
"file=\"" << output_filename <<
"\"/>" << std::endl;
9641 pvd_file <<
"</Collection>" << std::endl <<
"</VTKFile>";
////////////////////////////////////////////////////////////////////
void node_update(const bool &update_all_time_levels_for_new_node=false)
Broken assignment operator.
A class that contains the information required by Nodes that are located on Mesh boundaries....
std::map< unsigned, unsigned > *& index_of_first_value_assigned_by_face_element_pt()
Return pointer to the map giving the index of the first face element value.
GeneralisedTimestepper used to store the arclength derivatives and pervious solutions required in con...
void set_consistent_pinned_values(Data *const &data_pt)
Set consistent values of the derivatives and current value when the data is pinned....
void set_consistent_pinned_positions(Node *const &node_pt)
Set consistent values of the derivatives and current value when the Nodes position is pinned....
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
bool does_pointer_correspond_to_value(double *const ¶meter_pt)
Check whether the pointer parameter_pt addresses internal data values.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
bool is_halo() const
Is this Data a halo?
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
void dump(std::ostream &dump_file) const
Dump the data object to a file.
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
void set_nonhalo()
Label the node as not being a halo.
void set_halo(const unsigned &non_halo_proc_ID)
Label the node as halo and specify processor that holds non-halo counterpart.
void read(std::ifstream &restart_file)
Read data object from a file.
int non_halo_proc_ID()
ID of processor ID that holds non-halo counterpart of halo node; negative if not a halo.
Information for documentation of results: Directory and file number to enable output in the form RESL...
bool is_doc_enabled() const
Are we documenting?
std::string & label()
String used (e.g.) for labeling output files.
std::string directory() const
Output directory.
unsigned & number()
Number used (e.g.) for labeling output files.
A general Finite Element class.
virtual void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Broken virtual. Needs to be implemented for...
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output an exact solution over the element.
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
virtual std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
void check_jacobian(const double &jacobian) const
Helper function used to check for singular or negative Jacobians in the transform from local to globa...
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
unsigned nnode() const
Return the number of nodes.
virtual void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Return the paraview element type. Broken virtual. Needs to be implemented for each new geometric elem...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
void get_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
virtual unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Broken virtual. Needs to be implemented for each new...
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
virtual unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot. Broken virtual; can b...
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
virtual void scalar_value_fct_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
Write values of the i-th scalar field at the plot points. Broken virtual. Needs to be implemented for...
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
void output_paraview(std::ofstream &file_out, const unsigned &nplot) const
Paraview output – this outputs the coordinates at the plot points (for parameter nplot) to specified ...
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
static bool Suppress_output_while_checking_for_inverted_elements
Static boolean to suppress output while checking for inverted elements.
int get_node_number(Node *const &node_pt) const
Return the number of the node *node_pt if this node is in the element, else return -1;.
virtual void write_paraview_offsets(std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
Return the offsets for the paraview sub-elements. Broken virtual. Needs to be implemented for each ne...
virtual void write_paraview_output_offset_information(std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
Fill in the offset information for paraview plot. Broken virtual. Needs to be implemented for each ne...
A Generalised Element class.
bool is_halo() const
Is this element a halo?
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
unsigned ninternal_data() const
Return the number of internal data objects.
bool must_be_kept_as_halo() const
Test whether the element must be kept as a halo element.
void add_internal_value_pt_to_map(std::map< unsigned, double * > &map_of_value_pt)
Add pointers to the internal data values to map indexed by the global equation number.
Class that contains data for hanging nodes.
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
unsigned nmaster() const
Return the number of master nodes.
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Node * halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th halo node in this Mesh whose non-halo counterpart is held on processor p.
Vector< Node * > Node_pt
Vector of pointers to nodes.
std::map< unsigned, Vector< Node * > > Shared_node_pt
Map of vectors holding the pointers to the shared nodes. These are all the nodes that are on two "nei...
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
bool does_pointer_correspond_to_mesh_data(double *const ¶meter_pt)
Does the double pointer correspond to any mesh data.
unsigned ndof_types() const
Return number of dof types in mesh.
void remove_boundary_node(const unsigned &b, Node *const &node_pt)
Remove a node from the boundary b.
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
void output_external_halo_elements(std::ostream &outfile, const unsigned &n_plot=5)
Output all external halo elements.
void setup_shared_node_scheme()
Setup shared node scheme.
std::map< unsigned, Vector< GeneralisedElement * > > Root_haloed_element_pt
Map of vectors holding the pointers to the root haloed elements.
void add_halo_node_pt(const unsigned &p, Node *&nod_pt)
Add halo node whose non-halo counterpart is held on processor p to the storage scheme for halo nodes.
unsigned nexternal_haloed_element()
Total number of external haloed elements in this Mesh.
unsigned add_external_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Add external haloed node whose halo (external) counterpart is held on processor p to the storage sche...
Node * haloed_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th haloed node in this Mesh whose halo counterpart is held on processor p.
unsigned nexternal_halo_node()
Total number of external halo nodes in this Mesh.
unsigned nodal_dimension() const
Return number of nodal dimension in mesh.
bool Lookup_for_elements_next_boundary_is_setup
Flag to indicate that the lookup schemes for elements that are adjacent to the boundaries has been se...
void flush_element_and_node_storage()
Flush storage for elements and nodes by emptying the vectors that store the pointers to them....
std::map< unsigned, Vector< GeneralisedElement * > > External_haloed_element_pt
Map of vectors holding the pointers to the external haloed elements.
Node * shared_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th shared node in this Mesh who has a counterpart on processor p.
void add_shared_node_pt(const unsigned &p, Node *&nod_pt)
Add shared node whose counterpart is held on processor p to the storage scheme for shared nodes....
void flush_element_storage()
Flush storage for elements (only) by emptying the vectors that store the pointers to them....
Node *& external_halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external halo node in this Mesh whose non-halo external counterpart is held on...
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
void get_haloed_node_stats(double &av_number, unsigned &max_number, unsigned &min_number)
Get haloed node stats for this distributed mesh: Average/max/min number of haloed nodes over all proc...
void check_halo_schemes(DocInfo &doc_info, double &max_permitted_error_for_halo_check)
Check halo and shared schemes on the mesh.
virtual void set_mesh_level_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Function that can be used to set any additional timestepper data stored at the Mesh (as opposed to no...
Vector< GeneralisedElement * > haloed_element_pt(const unsigned &p)
Return vector of haloed elements in this Mesh whose haloing counterpart is held on processor p.
void describe_local_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the local dofs of the elements. The ostream specifies the output stream to which...
virtual void classify_halo_and_haloed_nodes(DocInfo &doc_info, const bool &report_stats)
Classify the halo and haloed nodes in the mesh. Virtual so it can be overloaded to perform additional...
void add_root_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add root haloed element whose non-halo counterpart is held on processor p to the storage scheme for h...
void synchronise_shared_nodes(const bool &report_stats)
Synchronise shared node lookup schemes to cater for the the case where: (1) a certain node on the cur...
void set_elemental_internal_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with the internal data stored within elements in the meah.
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
bool Keep_all_elements_as_halos
bool to indicate whether to keep all elements in a mesh as halos or not
std::map< unsigned, Vector< Node * > > External_haloed_node_pt
Map of vectors holding the pointers to the external haloed nodes.
virtual void setup_boundary_element_info()
Interface for function that is used to setup the boundary information (Empty virtual function – imple...
std::map< unsigned, Vector< GeneralisedElement * > > Root_halo_element_pt
Map of vectors holding the pointers to the root halo elements.
unsigned elemental_dimension() const
Return number of elemental dimension in mesh.
Node *& external_haloed_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external haloed node in this Mesh whose halo external counterpart is held on p...
void check_inverted_elements(bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
Check for inverted elements and report outcome in boolean variable. This visits all elements at their...
OomphCommunicator * Comm_pt
Pointer to communicator – set to NULL if mesh is not distributed.
void add_root_halo_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add root halo element whose non-halo counterpart is held on processor p to this Mesh.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt)
Output a given Vector function at f(n_plot) points in each element.
void shift_time_values()
Shift time-dependent data along for next timestep: Deal with nodal Data/positions and the element's i...
virtual void reorder_nodes(const bool &use_old_ordering=true)
Re-order nodes in the order in which they appear in elements – can be overloaded for more efficient r...
void prune_halo_elements_and_nodes(Vector< GeneralisedElement * > &deleted_element_pt, const bool &report_stats=false)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement,...
void calculate_predictions()
Calculate predictions for all Data and positions associated with the mesh, usually used in adaptive t...
unsigned nboundary() const
Return number of boundaries.
unsigned add_external_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add external haloed element whose non-halo counterpart is held on processor p to the storage scheme f...
void describe_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the dofs of the Mesh. The ostream specifies the output stream to which the descr...
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
void output_fct_paraview(std::ofstream &file_out, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
Output in paraview format into specified file. Breaks up each element into sub-elements for plotting ...
void remove_boundary_nodes()
Clear all pointers to boundary nodes.
bool Output_halo_elements
Bool for output of halo elements.
unsigned long nnode() const
Return number of nodes in the mesh.
virtual void node_update(const bool &update_all_solid_nodes=false)
Update nodal positions in response to changes in the domain shape. Uses the FiniteElement::get_x(....
void get_efficiency_of_mesh_distribution(double &av_efficiency, double &max_efficiency, double &min_efficiency)
Get efficiency of mesh distribution: In an ideal distribution without halo overhead,...
unsigned nroot_haloed_element()
Total number of root haloed elements in this Mesh.
void output_external_haloed_elements(std::ostream &outfile, const unsigned &n_plot=5)
Output all external haloed elements.
void delete_all_external_storage()
Wipe the storage for all externally-based elements.
void get_all_halo_data(std::map< unsigned, double * > &map_of_halo_data)
Get all the halo data stored in the mesh and add pointers to the data to the map, indexed by global e...
unsigned nhalo_node()
Total number of halo nodes in this Mesh.
std::map< unsigned, Vector< Node * > > Halo_node_pt
Map of vectors holding the pointers to the halo nodes.
virtual void create_shared_boundaries(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, const Vector< GeneralisedElement * > &backed_up_el_pt, const Vector< FiniteElement * > &backed_up_f_el_pt, std::map< Data *, std::set< unsigned >> &processors_associated_with_data, const bool &overrule_keep_as_halo_element_status)
Creates the shared boundaries, only used in unstructured meshes In this case with the "TriangleMesh" ...
unsigned check_for_repeated_nodes(const double &epsilon=1.0e-12)
Check for repeated nodes within a given spatial tolerance. Return (0/1) for (pass/fail).
void convert_to_boundary_node(Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
A function that upgrades an ordinary node to a boundary node We shouldn't ever really use this,...
void output(std::ostream &outfile)
Output for all elements.
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
std::map< unsigned, Vector< Node * > > Haloed_node_pt
Map of vectors holding the pointers to the haloed nodes.
void assign_initial_values_impulsive()
Assign initial values for an impulsive start.
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
bool Resize_halo_nodes_not_required
Set this to true to suppress resizing of halo nodes (at your own risk!)
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
void assign_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers in all elements If the boolean argument is true then also store poi...
Vector< Vector< Node * > > Boundary_node_pt
Vector of Vector of pointers to nodes on the boundaries: Boundary_node_pt(b,n). Note that this is pri...
virtual ~Mesh()
Virtual Destructor to clean up all memory.
void add_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Add haloed node whose halo counterpart is held on processor p to the storage scheme for haloed nodes.
void get_external_halo_node_pt(Vector< Node * > &external_halo_node_pt)
Get vector of pointers to all external halo nodes.
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
virtual void read(std::ifstream &restart_file)
Read solution from restart file.
void set_nodal_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with the nodal data in the mesh.
void doc_mesh_distribution(DocInfo &doc_info)
Doc the mesh distribution, to be processed with tecplot macros.
unsigned self_test()
Self-test: Check elements and nodes. Return 0 for OK.
void get_halo_node_stats(double &av_number, unsigned &max_number, unsigned &min_number)
Get halo node stats for this distributed mesh: Average/max/min number of halo nodes over all processo...
void set_consistent_pinned_values_for_continuation(ContinuationStorageScheme *const &continuation_stepper_pt)
Set consistent values for pinned data in continuation.
Vector< Node * > prune_dead_nodes()
Prune nodes. Nodes that have been marked as obsolete are removed from the mesh (and its boundary-node...
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
void resize_halo_nodes()
Helper function that resizes halo nodes to the same size as their non-halo counterparts if required....
std::map< unsigned, Vector< Node * > > External_halo_node_pt
Map of vectors holding the pointers to the external halo nodes.
unsigned nhaloed_node()
Total number of haloed nodes in this Mesh.
unsigned nexternal_haloed_node()
Total number of external haloed nodes in this Mesh.
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
Boolean used to control warning about empty mesh level timestepper function.
unsigned long assign_global_eqn_numbers(Vector< double * > &Dof_pt)
Assign the global equation numbers in the Data stored at the nodes and also internal element Data....
virtual void distribute(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, Vector< GeneralisedElement * > &deleted_element_pt, DocInfo &doc_info, const bool &report_stats, const bool &overrule_keep_as_halo_element_status)
Distribute the problem and doc; make this virtual to allow overloading for particular meshes where fu...
void null_external_halo_node(const unsigned &p, Node *nod_pt)
Null out specified external halo node (used when deleting duplicates)
unsigned nshared_node()
Total number of shared nodes in this Mesh.
void remove_null_pointers_from_external_halo_node_storage()
Consolidate external halo node storage by removing nulled out pointes in external halo and haloed sch...
unsigned long nelement() const
Return number of elements in the mesh.
virtual void get_node_reordering(Vector< Node * > &reordering, const bool &use_old_ordering=true) const
Get a reordering of the nodes in the order in which they appear in elements – can be overloaded for m...
void merge_meshes(const Vector< Mesh * > &sub_mesh_pt)
Merge meshes. Note: This simply merges the meshes' elements and nodes (ignoring duplicates; no bounda...
void output_paraview(std::ofstream &file_out, const unsigned &nplot) const
Output in paraview format into specified file. Breaks up each element into sub-elements for plotting ...
void output_boundaries(std::ostream &outfile)
Output the nodes on the boundaries (into separate tecplot zones)
unsigned nexternal_halo_element()
Total number of external halo elements in this Mesh.
virtual void dump(std::ofstream &dump_file, const bool &use_old_ordering=true) const
Dump the data in the mesh into a file for restart.
std::map< unsigned, Vector< GeneralisedElement * > > External_halo_element_pt
External halo(ed) elements are created as and when they are needed to act as source elements for the ...
Vector< GeneralisedElement * > halo_element_pt(const unsigned &p)
Return vector of halo elements in this Mesh whose non-halo counterpart is held on processor p.
unsigned nroot_halo_element()
Total number of root halo elements in this Mesh.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i). ‘Type’: k; Coordinate direction: i.
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
void copy(Node *orig_node_pt)
Copy all nodal data from specified Node object.
void set_nonhanging()
Label node as non-hanging node by removing all hanging node data.
bool is_obsolete()
Test whether node is obsolete.
void set_obsolete()
Mark node as obsolete.
void read(std::ifstream &restart_file)
Read nodal position and associated data from file for restart.
void position(Vector< double > &pos) const
Compute Vector of nodal positions either directly or via hanging node representation.
virtual void remove_from_boundary(const unsigned &b)
Broken interface for removing the node from the mesh boundary b Here to provide error reporting.
void set_non_obsolete()
Mark node as non-obsolete.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
virtual void add_to_boundary(const unsigned &b)
Broken interface for adding the node to the mesh boundary b Essentially here for error reporting.
unsigned hang_code()
Code that encapsulates the hanging status of the node (incl. the geometric hanging status) as .
void resize(const unsigned &n_value)
Resize the number of equations.
bool is_hanging() const
Test whether the node is geometrically hanging.
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
std::ostream *& stream_pt()
Access function for the stream pointer.
An OomphLibError object which should be thrown when an run-time error is encountered....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
An OomphLibWarning object which should be created as a temporary object to issue a warning....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
void get_father_at_refinement_level(unsigned &refinement_level, RefineableElement *&father_at_reflevel_pt)
Return a pointer to the "father" element at the specified refinement level.
unsigned refinement_level() const
Return the Refinement level.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void set_lagrangian_nodal_coordinates()
Make the current configuration the undeformed one by setting the nodal Lagrangian coordinates to thei...
static SolidICProblem Solid_IC_problem
Static problem that can be used to assign initial conditions on a given solid mesh (need to define th...
SolidNode * node_pt(const unsigned long &n)
Return a pointer to the n-th global SolidNode.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
unsigned nlagrangian_type() const
Number of types of Lagrangian coordinates used to interpolate the Lagrangian coordinates within the e...
double & xi_gen(const unsigned &k, const unsigned &i)
Reference to the generalised Lagrangian position. ‘Type’: k; 'Coordinate direction: i.
double lagrangian_position(const unsigned &i) const
Return lagrangian coordinate either directly or via hanging node representation.
void read(std::ifstream &restart_file)
Read nodal positions (variable and fixed) and associated data from file for restart.
double & xi(const unsigned &i)
Reference to i-th Lagrangian position.
void copy(SolidNode *orig_node_pt)
Copy nodal positions and associated data from specified node object.
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
unsigned nlagrangian() const
Return number of lagrangian coordinates.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
virtual void get_elements_at_refinement_level(unsigned &refinement_level, Vector< RefineableElement * > &level_elements)
Extract the elements at a particular refinement level in the refinement pattern - used in Mesh::redis...
unsigned uniform_refinement_level_when_pruned() const
Level to which the mesh was uniformly refined when it was pruned (const version)
virtual void get_refinement_levels(unsigned &min_refinement_level, unsigned &max_refinement_level)
Get max/min refinement levels in mesh.
virtual void setup_tree_forest()=0
Set up the tree forest associated with the Mesh (if any)
virtual void get_refinement_pattern(Vector< Vector< unsigned >> &to_be_refined)
Extract refinement pattern: Consider the hypothetical mesh obtained by truncating the refinement of t...
A generalised tree base class that abstracts the common functionality between the quad- and octrees u...
void stick_all_tree_nodes_into_vector(Vector< Tree * > &)
Traverse and stick pointers to all "nodes" into Vector.
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
bool is_leaf()
Return true if the tree is a leaf node.
void flush_object()
Flush the object represented by the tree.
void stick_leaves_into_vector(Vector< Tree * > &)
Traverse tree and stick pointers to leaf "nodes" (only) into Vector.
TreeRoot *& root_pt()
Return pointer to root of the tree.
Base class for triangle meshes (meshes made of 2D triangle elements). Note: we choose to template Tri...
virtual void reset_boundary_element_info(Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned >> &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
Virtual function to perform the reset boundary elements info rutines.
double region_attribute(const unsigned &i)
Return the attribute associated with region i.
unsigned nboundary_element_in_region(const unsigned &b, const unsigned &r) const
Return the number of elements adjacent to boundary b in region r.
unsigned nregion()
Return the number of regions specified by attributes.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
bool Doc_comprehensive_timings
Global boolean to switch on comprehensive timing – can probably be declared const false when developm...
bool node_global_position_comparison(Node *nd1_pt, Node *nd2_pt)
Function for ordering nodes. Return true if first node's position is "before" second nodes....
void write_pvd_footer(std::ofstream &pvd_file)
Write the pvd file footer.
void write_pvd_header(std::ofstream &pvd_file)
Write the pvd file header.
void write_pvd_information(std::ofstream &pvd_file, const std::string &output_filename, const double &time)
Add name of output file and associated continuous time to pvd file.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...