41 template<
unsigned INITIAL_NNODE_1D>
47 switch (this->nnode_1d())
74 oomph_info <<
"\n ERROR: Exceeded maximum polynomial order for";
75 oomph_info <<
"\n shape functions." << std::endl;
81 template<
unsigned INITIAL_NNODE_1D>
85 this->local_coordinate_of_node(n, s_fraction);
86 s_fraction[0] = 0.5 * (s_fraction[0] + 1.0);
89 template<
unsigned INITIAL_NNODE_1D>
91 const unsigned& n1d,
const unsigned&
i)
93 switch (this->nnode_1d())
120 std::ostringstream error_message;
121 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
122 error_message <<
"\n shape functions.\n";
124 OOMPH_CURRENT_FUNCTION,
125 OOMPH_EXCEPTION_LOCATION);
134 template<
unsigned INITIAL_NNODE_1D>
147 if (std::fabs(
s[0] + 1.0) < tol)
152 else if (std::fabs(
s[0] - 1.0) < tol)
154 index[0] = this->nnode_1d() - 1;
163 for (
unsigned n = 1; n < this->nnode_1d() - 1; n++)
165 if (std::fabs(z[n] -
s[0]) < tol)
175 return this->node_pt(index[0]);
185 template<
unsigned INITIAL_NNODE_1D>
200 template<
unsigned INITIAL_NNODE_1D>
202 Tree*
const& adopted_father_pt,
const unsigned& initial_p_order)
208 if (adopted_father_pt != 0)
211 father_pt =
dynamic_cast<BinaryTree*
>(adopted_father_pt);
214 else if (Tree_pt != 0)
217 father_pt =
dynamic_cast<BinaryTree*
>(binary_tree_pt()->father_pt());
227 if (father_pt != 0 || initial_p_order != 0)
234 if (father_el_pt != 0)
236 unsigned father_p_order = father_el_pt->
p_order();
238 P_order = father_p_order;
243 P_order = initial_p_order;
248 unsigned new_n_node = P_order;
251 this->set_n_node(new_n_node);
254 delete this->integral_pt();
276 std::ostringstream error_message;
277 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
278 error_message <<
"\n integration scheme.\n";
280 OOMPH_CURRENT_FUNCTION,
281 OOMPH_EXCEPTION_LOCATION);
290 template<
unsigned INITIAL_NNODE_1D>
326 template<
unsigned INITIAL_NNODE_1D>
335 if (clone_el_pt == 0)
338 "Cloned copy must be a PRefineableQElement<1,INITIAL_NNODE_1D>!",
339 OOMPH_CURRENT_FUNCTION,
340 OOMPH_EXCEPTION_LOCATION);
347 bool clone_is_ok =
true;
350 clone_is_ok = clone_is_ok && (clone_el_pt->
p_order() == this->p_order());
354 std::ostringstream err_stream;
355 err_stream <<
"Clone element has a different p-order from me,"
357 <<
"but it is supposed to be a copy!" << std::endl;
360 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
364 clone_is_ok = clone_is_ok && (clone_el_pt->
nnode() == this->nnode());
368 std::ostringstream err_stream;
369 err_stream <<
"Clone element has a different number of nodes from me,"
371 <<
"but it is supposed to be a copy!" << std::endl;
374 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
378 for (
unsigned n = 0; n < this->nnode(); n++)
381 clone_is_ok && (clone_el_pt->
node_pt(n) == this->node_pt(n));
386 std::ostringstream err_stream;
387 err_stream <<
"The nodes belonging to the clone element are different"
389 <<
"from mine, but it is supposed to be a copy!"
393 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
405 throw OomphLibError(
"Can't handle generalised nodal positions (yet).",
406 OOMPH_CURRENT_FUNCTION,
407 OOMPH_EXCEPTION_LOCATION);
412 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
415 unsigned ntstorage = time_stepper_pt->
ntstorage();
421 delete this->integral_pt();
443 std::ostringstream error_message;
444 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
445 error_message <<
"\n integration scheme.\n";
447 OOMPH_CURRENT_FUNCTION,
448 OOMPH_EXCEPTION_LOCATION);
452 this->set_n_node(P_order);
469 for (
unsigned n = 0; n < P_order; n++)
473 s_fraction[0] = this->local_one_d_fraction_of_node(n, 0);
476 s[0] = s_left[0] + (s_right[0] - s_left[0]) * s_fraction[0];
479 bool node_done =
false;
487 if (created_node_pt != 0)
490 this->node_pt(n) = created_node_pt;
498 for (
unsigned t = 0;
t < ntstorage;
t++)
507 unsigned n_val_at_node = created_node_pt->
nvalue();
508 unsigned n_val_from_function = prev_values.size();
510 unsigned n_var = n_val_at_node < n_val_from_function ?
514 for (
unsigned k = 0; k < n_var; k++)
516 created_node_pt->
set_value(
t, k, prev_values[k]);
535 created_node_pt = this->construct_node(n, time_stepper_pt);
550 for (
unsigned t = 0;
t < ntstorage;
t++)
557 clone_el_pt->
get_x(
t,
s, x_prev);
560 created_node_pt->
x(
t, 0) = x_prev[0];
574 const unsigned n_value = created_node_pt->
nvalue();
577 for (
unsigned v = 0; v < n_value; v++)
579 created_node_pt->
set_value(
t, v, prev_values[v]);
595 std::string error_message =
"Have not implemented p-refinement for";
596 error_message +=
"Algebraic p-refineable elements yet\n";
600 "PRefineableQElement<1,INITIAL_NNODE_1D>::p_refine()",
601 OOMPH_EXCEPTION_LOCATION);
625 this->node_pt(n),
s, clone_el_pt);
639 if (clone_m_el_pt != 0)
653 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
655 "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
656 error_message +=
"then I should be too....\n";
659 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
678 for (
unsigned n = 0; n < P_order; n++)
681 s_fraction[0] = this->local_one_d_fraction_of_node(n, 0);
683 s[0] = s_left[0] + (s_right[0] - s_left[0]) * s_fraction[0];
686 for (
unsigned t = 0;
t < ntstorage;
t++)
694 this->get_x(
t,
s, x_prev);
697 this->node_pt(n)->x(
t, 0) = x_prev[0];
705 this->further_build();
711 template<
unsigned INITIAL_NNODE_1D>
725 for (
unsigned i = 0;
i < p_order();
i++)
739 for (
unsigned i = 0;
i < p_order();
i++)
753 for (
unsigned i = 0;
i < p_order();
i++)
767 for (
unsigned i = 0;
i < p_order();
i++)
781 for (
unsigned i = 0;
i < p_order();
i++)
795 for (
unsigned i = 0;
i < p_order();
i++)
802 oomph_info <<
"\n ERROR: PRefineableQElement::shape() exceeded maximum";
803 oomph_info <<
"\n polynomial order for shape functions."
811 template<
unsigned INITIAL_NNODE_1D>
825 for (
unsigned i = 0;
i < p_order();
i++)
828 dpsi(
i, 0) = dpsi1ds[
i];
841 for (
unsigned i = 0;
i < p_order();
i++)
844 dpsi(
i, 0) = dpsi1ds[
i];
856 for (
unsigned i = 0;
i < p_order();
i++)
859 dpsi(
i, 0) = dpsi1ds[
i];
871 for (
unsigned i = 0;
i < p_order();
i++)
874 dpsi(
i, 0) = dpsi1ds[
i];
886 for (
unsigned i = 0;
i < p_order();
i++)
889 dpsi(
i, 0) = dpsi1ds[
i];
901 for (
unsigned i = 0;
i < p_order();
i++)
904 dpsi(
i, 0) = dpsi1ds[
i];
910 <<
"\n ERROR: PRefineableQElement::dshape_local() exceeded maximum";
911 oomph_info <<
"\n polynomial order for shape functions."
920 template<
unsigned INITIAL_NNODE_1D>
924 std::ostringstream error_message;
926 <<
"\nd2shape_local currently not implemented for this element\n";
928 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
935 template<
unsigned INITIAL_NNODE_1D>
937 const int& value_id,
const int& my_edge, std::ofstream& output_hangfile)
945 template<
unsigned INITIAL_NNODE_1D>
950 unsigned n_sons = this->tree_pt()->nsons();
952 unsigned max_son_p_order = 0;
953 for (
unsigned ison = 0; ison < n_sons; ison++)
956 this->tree_pt()->son_pt(ison)->object_pt());
957 son_p_order[ison] = el_pt->
p_order();
958 if (son_p_order[ison] > max_son_p_order)
959 max_son_p_order = son_p_order[ison];
962 unsigned old_Nnode = this->nnode();
963 unsigned old_P_order = this->p_order();
965 this->p_order() = max_son_p_order;
968 delete this->integral_pt();
969 switch (this->p_order())
990 std::ostringstream error_message;
991 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
992 error_message <<
"\n integration scheme.\n";
994 OOMPH_CURRENT_FUNCTION,
995 OOMPH_EXCEPTION_LOCATION);
1000 for (
unsigned n = 0; n < old_Nnode; n++)
1002 old_node_pt[n] = this->node_pt(n);
1006 this->set_n_node(this->p_order());
1012 this->node_pt(0) = old_node_pt[0];
1013 this->node_pt(this->p_order() - 1) = old_node_pt[old_P_order - 1];
1021 if (this->node_pt(0) == 0)
1024 OOMPH_CURRENT_FUNCTION,
1025 OOMPH_EXCEPTION_LOCATION);
1028 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
1031 const unsigned ntstorage = time_stepper_pt->
ntstorage();
1037 const unsigned n_node = this->nnode_1d();
1040 for (
unsigned n = 0; n < n_node; n++)
1043 s_fraction[0] = this->local_one_d_fraction_of_node(n, 0);
1046 s[0] = -1.0 + 2.0 * s_fraction[0];
1049 if (this->node_pt(n) == 0)
1052 bool is_periodic =
false;
1053 Node* created_node_pt =
1054 this->node_created_by_neighbour(s_fraction, is_periodic);
1057 if (created_node_pt != 0)
1063 OOMPH_CURRENT_FUNCTION,
1064 OOMPH_EXCEPTION_LOCATION);
1069 this->node_pt(n) = created_node_pt;
1079 using namespace BinaryTreeNames;
1082 if (s_fraction[0] < 0.5)
1085 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
1091 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
1098 this->tree_pt()->son_pt(son)->object_pt());
1104 if (n == 0 || n == n_node - 1)
1107 "I am trying to rebuild one of the (two) vertex nodes in\n";
1109 "this 1D element. It should not have been possible to delete\n";
1110 error_message +=
"either of these!\n";
1113 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1118 this->node_pt(n) = this->construct_node(n, time_stepper_pt);
1132 for (
unsigned t = 0;
t < ntstorage;
t++)
1139 son_el_pt->
get_x(
t, s_in_son, x_prev);
1142 this->node_pt(n)->x(
t, 0) = x_prev[0];
1156 const unsigned n_value = this->node_pt(n)->nvalue();
1159 for (
unsigned v = 0; v < n_value; v++)
1161 this->node_pt(n)->set_value(
t, v, prev_values[v]);
1182 "Have not implemented rebuilding from sons for";
1183 error_message +=
"Algebraic p-refineable elements yet\n";
1187 "PRefineableQElement<1,INITIAL_NNODE_1D>::rebuild_from_sons()",
1188 OOMPH_EXCEPTION_LOCATION);
1197 template<
unsigned INITIAL_NNODE_1D>
1209 template<
unsigned INITIAL_NNODE_1D>
1214 unsigned Nnode_1d = this->nnode_1d();
1215 unsigned n0 = n % Nnode_1d;
1216 unsigned n1 = n / Nnode_1d;
1251 std::ostringstream error_message;
1252 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1253 error_message <<
"\n shape functions.\n";
1255 OOMPH_CURRENT_FUNCTION,
1256 OOMPH_EXCEPTION_LOCATION);
1262 template<
unsigned INITIAL_NNODE_1D>
1266 this->local_coordinate_of_node(n, s_fraction);
1267 s_fraction[0] = 0.5 * (s_fraction[0] + 1.0);
1268 s_fraction[1] = 0.5 * (s_fraction[1] + 1.0);
1271 template<
unsigned INITIAL_NNODE_1D>
1273 const unsigned& n1d,
const unsigned&
i)
1275 switch (this->nnode_1d())
1302 std::ostringstream error_message;
1303 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1304 error_message <<
"\n shape functions.\n";
1306 OOMPH_CURRENT_FUNCTION,
1307 OOMPH_EXCEPTION_LOCATION);
1316 template<
unsigned INITIAL_NNODE_1D>
1320 unsigned Nnode_1d = this->nnode_1d();
1327 for (
unsigned i = 0;
i < 2;
i++)
1332 bool is_found =
false;
1335 if (std::fabs(
s[
i] + 1.0) < tol)
1341 else if (std::fabs(
s[
i] - 1.0) < tol)
1343 index[
i] = Nnode_1d - 1;
1353 for (
unsigned n = 1; n < Nnode_1d - 1; n++)
1355 if (std::fabs(z[n] -
s[
i]) < tol)
1371 return this->node_pt(index[0] + Nnode_1d * index[1]);
1381 template<
unsigned INITIAL_NNODE_1D>
1385 using namespace QuadTreeNames;
1389 if (s_fraction[0] == 0.0)
1393 if (s_fraction[0] == 1.0)
1397 if (s_fraction[1] == 0.0)
1401 if (s_fraction[1] == 1.0)
1407 unsigned n_size = edges.size();
1419 int neigh_edge, diff_level;
1420 bool in_neighbouring_tree;
1423 for (
unsigned j = 0; j < n_size; j++)
1433 in_neighbouring_tree);
1442 bool a_vertex_node_is_built =
false;
1445 if (neigh_obj_pt == 0)
1448 "PRefineableQElement<2,INITIAL_NNODE_1D>::node_"
1449 "created_by_neighbour()",
1450 OOMPH_EXCEPTION_LOCATION);
1452 for (
unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node(); vnode++)
1454 if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
1455 a_vertex_node_is_built =
true;
1458 if (a_vertex_node_is_built)
1467 for (
unsigned i = 0;
i < 2;
i++)
1469 s[
i] = s_lo_neigh[
i] +
1470 s_fraction[translate_s[
i]] * (s_hi_neigh[
i] - s_lo_neigh[
i]);
1474 Node* neighbour_node_pt =
1478 if (neighbour_node_pt != 0)
1482 if (in_neighbouring_tree)
1486 quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1489 return neighbour_node_pt;
1505 template<
unsigned INITIAL_NNODE_1D>
1510 using namespace QuadTreeNames;
1514 if (s_fraction[0] == 0.0)
1518 if (s_fraction[0] == 1.0)
1522 if (s_fraction[1] == 0.0)
1526 if (s_fraction[1] == 1.0)
1532 unsigned n_size = edges.size();
1544 int neigh_edge, diff_level;
1545 bool in_neighbouring_tree;
1548 for (
unsigned j = 0; j < n_size; j++)
1558 in_neighbouring_tree);
1574 for (
unsigned i = 0;
i < 2;
i++)
1576 s[
i] = s_lo_neigh[
i] +
1577 s_fraction[translate_s[
i]] * (s_hi_neigh[
i] - s_lo_neigh[
i]);
1581 if (neigh_pt->
nsons() != 0)
1596 s_in_son[0] = 1.0 + 2.0 *
s[0];
1597 s_in_son[1] = 1.0 + 2.0 *
s[1];
1604 s_in_son[0] = 1.0 + 2.0 *
s[0];
1605 s_in_son[1] = -1.0 + 2.0 *
s[1];
1615 s_in_son[0] = -1.0 + 2.0 *
s[0];
1616 s_in_son[1] = 1.0 + 2.0 *
s[1];
1623 s_in_son[0] = -1.0 + 2.0 *
s[0];
1624 s_in_son[1] = -1.0 + 2.0 *
s[1];
1629 Node* neighbour_son_node_pt =
1634 if (neighbour_son_node_pt != 0)
1638 if (in_neighbouring_tree)
1642 quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1645 return neighbour_son_node_pt;
1666 template<
unsigned INITIAL_NNODE_1D>
1668 Tree*
const& adopted_father_pt,
const unsigned& initial_p_order)
1674 if (adopted_father_pt != 0)
1677 father_pt =
dynamic_cast<QuadTree*
>(adopted_father_pt);
1680 else if (Tree_pt != 0)
1683 father_pt =
dynamic_cast<QuadTree*
>(quadtree_pt()->father_pt());
1693 if (father_pt != 0 || initial_p_order != 0)
1700 if (father_el_pt != 0)
1702 unsigned father_p_order = father_el_pt->
p_order();
1704 P_order = father_p_order;
1709 P_order = initial_p_order;
1714 unsigned new_n_node = P_order * P_order;
1717 this->set_n_node(new_n_node);
1720 delete this->integral_pt();
1742 std::ostringstream error_message;
1743 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1744 error_message <<
"\n integration scheme.\n";
1746 OOMPH_CURRENT_FUNCTION,
1747 OOMPH_EXCEPTION_LOCATION);
1756 template<
unsigned INITIAL_NNODE_1D>
1760 using namespace QuadTreeNames;
1763 unsigned n_p = nnode_1d();
1766 if (Father_bound[n_p].nrow() == 0)
1768 setup_father_bounds();
1775 int son_type = Tree_pt->
son_type();
1784 "Something fishy here: I have no father and yet \n";
1785 error_message +=
"I have no nodes. Who has created me then?!\n";
1788 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1802 unsigned ntstorage = time_stepper_pt->ntstorage();
1848 for (
unsigned i = 0;
i < 2;
i++)
1852 0.5 * (s_lo[
i] + 1.0) *
1856 0.5 * (s_hi[
i] + 1.0) *
1863 if (father_el_pt->
node_pt(0) == 0)
1866 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1867 OOMPH_CURRENT_FUNCTION,
1868 OOMPH_EXCEPTION_LOCATION);
1878 for (
unsigned i0 = 0; i0 < n_p; i0++)
1881 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
1883 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
1885 for (
unsigned i1 = 0; i1 < n_p; i1++)
1888 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
1890 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
1893 jnod = i0 + n_p * i1;
1914 Node* created_node_pt =
1919 if (created_node_pt != 0)
1922 this->node_pt(jnod) = created_node_pt;
1930 for (
unsigned t = 0;
t < ntstorage;
t++)
1939 unsigned n_val_at_node = created_node_pt->
nvalue();
1940 unsigned n_val_from_function = prev_values.size();
1942 unsigned n_var = n_val_at_node < n_val_from_function ?
1944 n_val_from_function;
1946 for (
unsigned k = 0; k < n_var; k++)
1948 created_node_pt->
set_value(
t, k, prev_values[k]);
1989 template<
unsigned INITIAL_NNODE_1D>
2007 unsigned ngeom_object = m_el_pt->ngeom_object();
2009 for (
unsigned i = 0;
i < ngeom_object;
i++)
2018 if (clone_el_pt == 0)
2021 "Cloned copy must be a PRefineableQElement<2,INITIAL_NNODE_1D>!",
2022 OOMPH_CURRENT_FUNCTION,
2023 OOMPH_EXCEPTION_LOCATION);
2030 bool clone_is_ok =
true;
2033 clone_is_ok = clone_is_ok && (clone_el_pt->
p_order() == this->p_order());
2037 std::ostringstream err_stream;
2038 err_stream <<
"Clone element has a different p-order from me,"
2040 <<
"but it is supposed to be a copy!" << std::endl;
2043 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2047 clone_is_ok = clone_is_ok && (clone_el_pt->
nnode() == this->nnode());
2051 std::ostringstream err_stream;
2052 err_stream <<
"Clone element has a different number of nodes from me,"
2054 <<
"but it is supposed to be a copy!" << std::endl;
2057 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2061 for (
unsigned n = 0; n < this->nnode(); n++)
2064 clone_is_ok && (clone_el_pt->
node_pt(n) == this->node_pt(n));
2069 std::ostringstream err_stream;
2070 err_stream <<
"The nodes belonging to the clone element are different"
2072 <<
"from mine, but it is supposed to be a copy!"
2076 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2094 clone_is_ok && (geom_object_pt.size() == clone_geom_object_pt.size());
2095 for (
unsigned n = 0;
2096 n < std::min(geom_object_pt.size(), clone_geom_object_pt.size());
2100 clone_is_ok && (clone_geom_object_pt[n] == geom_object_pt[n]);
2105 std::ostringstream err_stream;
2106 err_stream <<
"The clone element has different geometric objects"
2108 <<
"from mine, but it is supposed to be a copy!"
2113 "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2114 OOMPH_EXCEPTION_LOCATION);
2127 throw OomphLibError(
"Can't handle generalised nodal positions (yet).",
2128 OOMPH_CURRENT_FUNCTION,
2129 OOMPH_EXCEPTION_LOCATION);
2134 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
2137 unsigned ntstorage = time_stepper_pt->
ntstorage();
2143 delete this->integral_pt();
2165 std::ostringstream error_message;
2166 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
2167 error_message <<
"\n integration scheme.\n";
2169 OOMPH_CURRENT_FUNCTION,
2170 OOMPH_EXCEPTION_LOCATION);
2174 this->set_n_node(P_order * P_order);
2195 for (
unsigned i0 = 0; i0 < P_order; i0++)
2198 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
2200 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
2202 for (
unsigned i1 = 0; i1 < P_order; i1++)
2205 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
2207 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
2210 jnod = i0 + P_order * i1;
2214 bool node_done =
false;
2222 if (created_node_pt != 0)
2225 this->node_pt(jnod) = created_node_pt;
2233 for (
unsigned t = 0;
t < ntstorage;
t++)
2242 unsigned n_val_at_node = created_node_pt->
nvalue();
2243 unsigned n_val_from_function = prev_values.size();
2245 unsigned n_var = n_val_at_node < n_val_from_function ?
2247 n_val_from_function;
2249 for (
unsigned k = 0; k < n_var; k++)
2251 created_node_pt->
set_value(
t, k, prev_values[k]);
2267 bool is_periodic =
false;
2268 created_node_pt = node_created_by_neighbour(s_fraction, is_periodic);
2271 if (created_node_pt != 0)
2280 Node* neighbour_node_pt = created_node_pt;
2286 else if (s_fraction[0] == 1.0)
2288 else if (s_fraction[1] == 0.0)
2290 else if (s_fraction[1] == 1.0)
2297 std::set<unsigned> boundaries;
2309 if (boundaries.size() > 1)
2312 "boundaries.size()!=1 seems a bit strange..\n",
2313 OOMPH_CURRENT_FUNCTION,
2314 OOMPH_EXCEPTION_LOCATION);
2318 if (boundaries.size() == 0)
2320 std::ostringstream error_stream;
2321 error_stream <<
"Periodic node is not on a boundary...\n"
2322 <<
"Coordinates: " << created_node_pt->
x(0) <<
" "
2323 << created_node_pt->
x(1) <<
"\n";
2325 OOMPH_CURRENT_FUNCTION,
2326 OOMPH_EXCEPTION_LOCATION);
2332 this->construct_boundary_node(jnod, time_stepper_pt);
2337 for (
unsigned t = 0;
t < ntstorage;
t++)
2345 clone_el_pt->
get_x(
t,
s, x_prev);
2347 for (
unsigned i = 0;
i < 2;
i++)
2349 created_node_pt->
x(
t,
i) = x_prev[
i];
2358 for (std::set<unsigned>::iterator it = boundaries.begin();
2359 it != boundaries.end();
2373 *it, my_bound,
s, zeta);
2387 this->node_pt(jnod) = created_node_pt;
2401 bool is_periodic =
false;
2403 node_created_by_son_of_neighbour(s_fraction, is_periodic);
2406 if (created_node_pt != 0)
2415 Node* neighbour_node_pt = created_node_pt;
2421 else if (s_fraction[0] == 1.0)
2423 else if (s_fraction[1] == 0.0)
2425 else if (s_fraction[1] == 1.0)
2432 std::set<unsigned> boundaries;
2444 if (boundaries.size() > 1)
2447 "boundaries.size()!=1 seems a bit strange..\n",
2448 OOMPH_CURRENT_FUNCTION,
2449 OOMPH_EXCEPTION_LOCATION);
2453 if (boundaries.size() == 0)
2455 std::ostringstream error_stream;
2456 error_stream <<
"Periodic node is not on a boundary...\n"
2457 <<
"Coordinates: " << created_node_pt->
x(0)
2458 <<
" " << created_node_pt->
x(1) <<
"\n";
2460 OOMPH_CURRENT_FUNCTION,
2461 OOMPH_EXCEPTION_LOCATION);
2467 this->construct_boundary_node(jnod, time_stepper_pt);
2472 for (
unsigned t = 0;
t < ntstorage;
t++)
2480 clone_el_pt->
get_x(
t,
s, x_prev);
2482 for (
unsigned i = 0;
i < 2;
i++)
2484 created_node_pt->
x(
t,
i) = x_prev[
i];
2493 for (std::set<unsigned>::iterator it = boundaries.begin();
2494 it != boundaries.end();
2508 *it, my_bound,
s, zeta);
2522 this->node_pt(jnod) = created_node_pt;
2541 else if (s_fraction[0] == 1.0)
2543 else if (s_fraction[1] == 0.0)
2545 else if (s_fraction[1] == 1.0)
2552 std::set<unsigned> boundaries;
2564 if (boundaries.size() > 1)
2566 throw OomphLibError(
"boundaries.size()!=1 seems a bit strange..\n",
2567 OOMPH_CURRENT_FUNCTION,
2568 OOMPH_EXCEPTION_LOCATION);
2574 if (boundaries.size() > 0)
2578 this->construct_boundary_node(jnod, time_stepper_pt);
2585 Vector<int> bound_cons(this->ncont_interpolated_values());
2586 clone_el_pt->
get_bcs(my_bound, bound_cons);
2589 unsigned n_value = created_node_pt->
nvalue();
2590 for (
unsigned k = 0; k < n_value; k++)
2594 created_node_pt->
pin(k);
2601 dynamic_cast<SolidNode*
>(created_node_pt);
2602 if (solid_node_pt != 0)
2605 unsigned n_dim = created_node_pt->
ndim();
2610 if (clone_solid_el_pt == 0)
2613 "We have a SolidNode outside a refineable SolidElement\n";
2615 "during mesh refinement -- this doesn't make sense";
2618 OOMPH_CURRENT_FUNCTION,
2619 OOMPH_EXCEPTION_LOCATION);
2622 clone_solid_el_pt->
get_solid_bcs(my_bound, solid_bound_cons);
2625 for (
unsigned k = 0; k < n_dim; k++)
2627 if (solid_bound_cons[k])
2640 for (std::set<unsigned>::iterator it = boundaries.begin();
2641 it != boundaries.end();
2655 *it, my_bound,
s, zeta);
2667 created_node_pt = this->construct_node(jnod, time_stepper_pt);
2683 for (
unsigned t = 0;
t < ntstorage;
t++)
2691 clone_el_pt->
get_x(
t,
s, x_prev);
2694 for (
unsigned i = 0;
i < 2;
i++)
2696 created_node_pt->
x(
t,
i) = x_prev[
i];
2701 for (
unsigned t = 0;
t < ntstorage;
t++)
2708 unsigned n_value = created_node_pt->
nvalue();
2709 for (
unsigned k = 0; k < n_value; k++)
2711 created_node_pt->
set_value(
t, k, prev_values[k]);
2727 std::string error_message =
"Have not implemented p-refinement for";
2728 error_message +=
"Algebraic p-refineable elements yet\n";
2732 "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2733 OOMPH_EXCEPTION_LOCATION);
2757 this->node_pt(jnod),
s, clone_el_pt);
2773 if (clone_m_el_pt != 0)
2787 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
2789 "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
2790 error_message +=
"then I should be too....\n";
2793 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2858 for (
unsigned i0 = 0; i0 < P_order; i0++)
2861 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
2863 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
2865 for (
unsigned i1 = 0; i1 < P_order; i1++)
2868 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
2870 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
2873 jnod = i0 + P_order * i1;
2876 for (
unsigned t = 0;
t < ntstorage;
t++)
2884 this->get_x(
t,
s, x_prev);
2887 for (
unsigned i = 0;
i < 2;
i++)
2889 this->node_pt(jnod)->x(
t,
i) = x_prev[
i];
2900 this->further_build();
2906 template<
unsigned INITIAL_NNODE_1D>
2920 for (
unsigned i = 0;
i < 2;
i++)
2922 for (
unsigned j = 0; j < 2; j++)
2924 psi(2 *
i + j) = psi2[
i] * psi1[j];
2937 for (
unsigned i = 0;
i < 3;
i++)
2939 for (
unsigned j = 0; j < 3; j++)
2941 psi(3 *
i + j) = psi2[
i] * psi1[j];
2954 for (
unsigned i = 0;
i < 4;
i++)
2956 for (
unsigned j = 0; j < 4; j++)
2958 psi(4 *
i + j) = psi2[
i] * psi1[j];
2971 for (
unsigned i = 0;
i < 5;
i++)
2973 for (
unsigned j = 0; j < 5; j++)
2975 psi(5 *
i + j) = psi2[
i] * psi1[j];
2988 for (
unsigned i = 0;
i < 6;
i++)
2990 for (
unsigned j = 0; j < 6; j++)
2992 psi(6 *
i + j) = psi2[
i] * psi1[j];
3005 for (
unsigned i = 0;
i < 7;
i++)
3007 for (
unsigned j = 0; j < 7; j++)
3009 psi(7 *
i + j) = psi2[
i] * psi1[j];
3015 std::ostringstream error_message;
3016 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
3017 error_message <<
"\n polynomial order for shape functions.\n";
3019 OOMPH_CURRENT_FUNCTION,
3020 OOMPH_EXCEPTION_LOCATION);
3027 template<
unsigned INITIAL_NNODE_1D>
3043 for (
unsigned i = 0;
i < 2;
i++)
3045 for (
unsigned j = 0; j < 2; j++)
3048 dpsids(index, 0) = psi2[
i] * dpsi1ds[j];
3049 dpsids(index, 1) = dpsi2ds[
i] * psi1[j];
3050 psi[index] = psi2[
i] * psi1[j];
3067 for (
unsigned i = 0;
i < 3;
i++)
3069 for (
unsigned j = 0; j < 3; j++)
3072 dpsids(index, 0) = psi2[
i] * dpsi1ds[j];
3073 dpsids(index, 1) = dpsi2ds[
i] * psi1[j];
3074 psi[index] = psi2[
i] * psi1[j];
3091 for (
unsigned i = 0;
i < 4;
i++)
3093 for (
unsigned j = 0; j < 4; j++)
3096 dpsids(index, 0) = psi2[
i] * dpsi1ds[j];
3097 dpsids(index, 1) = dpsi2ds[
i] * psi1[j];
3098 psi[index] = psi2[
i] * psi1[j];
3115 for (
unsigned i = 0;
i < 5;
i++)
3117 for (
unsigned j = 0; j < 5; j++)
3120 dpsids(index, 0) = psi2[
i] * dpsi1ds[j];
3121 dpsids(index, 1) = dpsi2ds[
i] * psi1[j];
3122 psi[index] = psi2[
i] * psi1[j];
3139 for (
unsigned i = 0;
i < 6;
i++)
3141 for (
unsigned j = 0; j < 6; j++)
3144 dpsids(index, 0) = psi2[
i] * dpsi1ds[j];
3145 dpsids(index, 1) = dpsi2ds[
i] * psi1[j];
3146 psi[index] = psi2[
i] * psi1[j];
3163 for (
unsigned i = 0;
i < 7;
i++)
3165 for (
unsigned j = 0; j < 7; j++)
3168 dpsids(index, 0) = psi2[
i] * dpsi1ds[j];
3169 dpsids(index, 1) = dpsi2ds[
i] * psi1[j];
3170 psi[index] = psi2[
i] * psi1[j];
3178 std::ostringstream error_message;
3179 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
3180 error_message <<
"\n polynomial order for shape functions.\n";
3182 OOMPH_CURRENT_FUNCTION,
3183 OOMPH_EXCEPTION_LOCATION);
3191 template<
unsigned INITIAL_NNODE_1D>
3195 std::ostringstream error_message;
3197 <<
"\nd2shape_local currently not implemented for this element\n";
3199 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3206 template<
unsigned INITIAL_NNODE_1D>
3210 using namespace QuadTreeNames;
3213 unsigned n_sons = this->tree_pt()->nsons();
3215 unsigned max_son_p_order = 0;
3216 for (
unsigned ison = 0; ison < n_sons; ison++)
3219 this->tree_pt()->son_pt(ison)->object_pt());
3220 son_p_order[ison] = el_pt->
p_order();
3221 if (son_p_order[ison] > max_son_p_order)
3222 max_son_p_order = son_p_order[ison];
3225 unsigned old_Nnode = this->nnode();
3226 unsigned old_P_order = this->p_order();
3228 this->p_order() = max_son_p_order;
3231 delete this->integral_pt();
3232 switch (this->p_order())
3253 std::ostringstream error_message;
3254 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
3255 error_message <<
"\n integration scheme.\n";
3257 OOMPH_CURRENT_FUNCTION,
3258 OOMPH_EXCEPTION_LOCATION);
3263 for (
unsigned n = 0; n < old_Nnode; n++)
3265 old_node_pt[n] = this->node_pt(n);
3269 this->set_n_node(this->p_order() * this->p_order());
3272 this->node_pt(0) = old_node_pt[0];
3273 this->node_pt(this->p_order() - 1) = old_node_pt[old_P_order - 1];
3274 this->node_pt(this->p_order() * (this->p_order() - 1)) =
3275 old_node_pt[(old_P_order) * (old_P_order - 1)];
3276 this->node_pt(this->p_order() * this->p_order() - 1) =
3277 old_node_pt[(old_P_order) * (old_P_order)-1];
3280 if (this->p_order() % 2 == 1)
3283 unsigned midpoint = (this->p_order() - 1) / 2;
3287 quadtree_pt()->son_pt(
SW)->object_pt())
3288 ->vertex_node_pt(1);
3290 this->node_pt(midpoint * this->p_order()) =
3292 quadtree_pt()->son_pt(
SW)->object_pt())
3293 ->vertex_node_pt(2);
3295 this->node_pt((this->p_order() - 1) * this->p_order() + midpoint) =
3297 quadtree_pt()->son_pt(
NE)->object_pt())
3298 ->vertex_node_pt(2);
3300 this->node_pt((midpoint + 1) * this->p_order() - 1) =
3302 quadtree_pt()->son_pt(
NE)->object_pt())
3303 ->vertex_node_pt(1);
3309 if (this->node_pt(0) == 0)
3312 OOMPH_CURRENT_FUNCTION,
3313 OOMPH_EXCEPTION_LOCATION);
3316 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
3317 unsigned ntstorage = time_stepper_pt->
ntstorage();
3322 unsigned n_p = this->nnode_1d();
3323 for (
unsigned i0 = 0; i0 < n_p; i0++)
3326 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3328 s[0] = -1.0 + 2.0 * s_fraction[0];
3330 for (
unsigned i1 = 0; i1 < n_p; i1++)
3333 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
3335 s[1] = -1.0 + 2.0 * s_fraction[1];
3338 jnod = i0 + n_p * i1;
3342 bool node_done =
false;
3346 Node* created_node_pt = this->get_node_at_local_coordinate(
s);
3350 if (created_node_pt != 0)
3353 this->node_pt(jnod) = created_node_pt;
3367 bool is_periodic =
false;
3368 created_node_pt = node_created_by_neighbour(s_fraction, is_periodic);
3371 if (created_node_pt != 0)
3377 OOMPH_CURRENT_FUNCTION,
3378 OOMPH_EXCEPTION_LOCATION);
3383 this->node_pt(jnod) = created_node_pt;
3397 using namespace QuadTreeNames;
3400 if (s_fraction[0] < 0.5)
3403 if (s_fraction[1] < 0.5)
3407 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
3408 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
3415 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
3416 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
3422 if (s_fraction[1] < 0.5)
3426 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
3427 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
3434 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
3435 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
3443 this->tree_pt()->son_pt(son)->object_pt());
3455 else if (i0 == n_p - 1)
3479 else if (i1 == n_p - 1)
3497 std::set<unsigned> boundaries;
3511 if (boundaries.size() > 0)
3515 this->construct_boundary_node(jnod, time_stepper_pt);
3518 Vector<int> bound_cons(this->ncont_interpolated_values());
3519 son_el_pt->
get_bcs(boundary, bound_cons);
3522 unsigned nval = created_node_pt->
nvalue();
3523 for (
unsigned k = 0; k < nval; k++)
3527 created_node_pt->
pin(k);
3534 dynamic_cast<SolidNode*
>(created_node_pt);
3535 if (solid_node_pt != 0)
3538 unsigned n_dim = created_node_pt->
ndim();
3543 if (son_solid_el_pt == 0)
3546 "We have a SolidNode outside a refineable SolidElement\n";
3548 "during mesh refinement -- this doesn't make sense\n";
3551 OOMPH_CURRENT_FUNCTION,
3552 OOMPH_EXCEPTION_LOCATION);
3558 for (
unsigned k = 0; k < n_dim; k++)
3560 if (solid_bound_cons[k])
3570 for (std::set<unsigned>::iterator it = boundaries.begin();
3571 it != boundaries.end();
3585 *it, boundary, s_in_son, zeta);
3596 created_node_pt = this->construct_node(jnod, time_stepper_pt);
3613 for (
unsigned t = 0;
t < ntstorage;
t++)
3615 using namespace QuadTreeNames;
3620 son_el_pt->
get_x(
t, s_in_son, x_prev);
3621 for (
unsigned i = 0;
i < 2;
i++)
3623 created_node_pt->
x(
t,
i) = x_prev[
i];
3629 for (
unsigned t = 0;
t < ntstorage;
t++)
3637 for (
unsigned k = 0; k < created_node_pt->
nvalue(); k++)
3639 created_node_pt->
set_value(
t, k, prev_values[k]);
3654 "Have not implemented rebuilding from sons for";
3655 error_message +=
"Algebraic p-refineable elements yet\n";
3659 "PRefineableQElement<2,INITIAL_NNODE_1D>::rebuild_from_sons()",
3660 OOMPH_EXCEPTION_LOCATION);
3677 for (
unsigned j = 0; j < this->nnode(); j++)
3681 this->local_coordinate_of_node(j, s_in_node_update_element);
3688 ->set_node_update_info(
3689 this, s_in_node_update_element, geom_object_pt);
3701 for (
unsigned i0 = 0; i0 < n_p; i0++)
3704 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3706 s[0] = -1.0 + 2.0 * s_fraction[0];
3708 for (
unsigned i1 = 0; i1 < n_p; i1++)
3711 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
3713 s[1] = -1.0 + 2.0 * s_fraction[1];
3716 jnod = i0 + n_p * i1;
3719 for (
unsigned t = 0;
t < ntstorage;
t++)
3727 this->get_x(
t,
s, x_prev);
3730 for (
unsigned i = 0;
i < 2;
i++)
3732 this->node_pt(jnod)->x(
t,
i) = x_prev[
i];
3747 template<
unsigned INITIAL_NNODE_1D>
3761 if (this->macro_elem_pt() != 0)
3767 using namespace QuadTreeNames;
3770 unsigned n_p = nnode_1d();
3774 unsigned n_time = 1;
3779 double max_error_val = 0.0;
3788 for (
unsigned edge_counter = 0; edge_counter < 4; edge_counter++)
3792 int neigh_edge, diff_level;
3793 bool in_neighbouring_tree;
3803 in_neighbouring_tree);
3810 bool is_periodic =
false;
3811 if (in_neighbouring_tree)
3815 edges[edge_counter]);
3821 bool exclude_this_edge =
false;
3822 if (diff_level != 0)
3825 exclude_this_edge =
true;
3827 else if (neigh_pt->
nsons() != 0)
3830 exclude_this_edge =
true;
3834 unsigned my_p_order = this->p_order();
3835 unsigned neigh_p_order =
3838 if (my_p_order != neigh_p_order)
3841 exclude_this_edge =
true;
3856 if (!exclude_this_edge)
3859 for (
unsigned i0 = 0; i0 < n_p; i0++)
3862 Node* local_node_pt = 0;
3864 switch (edge_counter)
3868 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3869 s_fraction[1] = 0.0;
3871 local_node_pt = this->node_pt(i0);
3876 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3877 s_fraction[1] = 1.0;
3879 local_node_pt = this->node_pt(i0 + n_p * (n_p - 1));
3884 s_fraction[0] = 0.0;
3885 s_fraction[1] = this->local_one_d_fraction_of_node(i0, 1);
3887 local_node_pt = this->node_pt(n_p * i0);
3892 s_fraction[0] = 1.0;
3893 s_fraction[1] = this->local_one_d_fraction_of_node(i0, 1);
3895 local_node_pt = this->node_pt(n_p - 1 + n_p * i0);
3902 for (
unsigned i = 0;
i < 2;
i++)
3905 s[
i] = -1.0 + 2.0 * s_fraction[
i];
3909 s_fraction[translate_s[
i]] * (s_hi_neigh[
i] - s_lo_neigh[
i]);
3913 for (
unsigned t = 0;
t < n_time;
t++)
3918 t, s_in_neighb, x_in_neighb);
3921 if (is_periodic ==
false)
3923 for (
int i = 0;
i < 2;
i++)
3927 std::fabs(local_node_pt->
x(
t,
i) - x_in_neighb[
i]);
3933 << local_node_pt->
x(
t,
i) <<
" "
3934 << x_in_neighb[
i] << std::endl;
3937 << local_node_pt->
x(1) << std::endl;
3942 if (err > max_error_x[
i])
3944 max_error_x[
i] = err;
3954 t, s_in_neighb, values_in_neighb);
3958 this->get_interpolated_values(
t,
s, values);
3966 for (
unsigned ival = 0; ival < num_val; ival++)
3968 double err = std::fabs(values[ival] - values_in_neighb[ival]);
3973 << local_node_pt->
x(1) <<
" \n# "
3974 <<
"erru (S)" << err <<
" " << ival <<
" "
3975 << this->get_node_number(local_node_pt) <<
" "
3976 << values[ival] <<
" " << values_in_neighb[ival]
3980 if (err > max_error_val)
3982 max_error_val = err;
3991 max_error = max_error_x[0];
3992 if (max_error_x[1] > max_error) max_error = max_error_x[1];
3993 if (max_error_val > max_error) max_error = max_error_val;
3995 if (max_error > 1
e-9)
3997 oomph_info <<
"\n#------------------------------------ \n#Max error ";
3998 oomph_info << max_error_x[0] <<
" " << max_error_x[1] <<
" "
3999 << max_error_val << std::endl;
4000 oomph_info <<
"#------------------------------------ \n " << std::endl;
4050 template<
unsigned INITIAL_NNODE_1D>
4052 const int& value_id,
const int& my_edge, std::ofstream& output_hangfile)
4054 using namespace QuadTreeNames;
4059 int neigh_edge, diff_level;
4060 bool in_neighbouring_tree;
4070 in_neighbouring_tree);
4077 bool h_type_dependent =
false;
4079 bool p_type_dependent =
false;
4085 if (diff_level != 0)
4088 h_type_dependent =
true;
4091 else if (neigh_pt->
nsons() == 0)
4095 unsigned my_p_order =
4098 unsigned neigh_p_order =
4104 if (neigh_p_order == my_p_order)
4108 else if (neigh_p_order < my_p_order)
4111 p_type_dependent =
true;
4133 if (h_type_dependent || p_type_dependent)
4136 unsigned active_coord_index;
4137 if (my_edge ==
N || my_edge ==
S) active_coord_index = 0;
4138 else if (my_edge ==
E || my_edge ==
W)
4139 active_coord_index = 1;
4143 OOMPH_CURRENT_FUNCTION,
4144 OOMPH_EXCEPTION_LOCATION);
4159 const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
4164 bool is_periodic =
false;
4165 if (in_neighbouring_tree)
4168 this->tree_pt()->root_pt()->is_neighbour_periodic(my_edge);
4180 "Cannot do mortaring with periodic hanging nodes yet! (Fix me!)",
4181 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4182 OOMPH_EXCEPTION_LOCATION);
4187 unsigned neighbour_node_number = 0;
4188 Node* neighbour_node_pt = 0;
4191 for (
unsigned i0 = 0; i0 < neigh_n_p; i0++)
4197 neighbour_node_number = i0 + neigh_n_p * (neigh_n_p - 1);
4199 neighbour_node_number, value_id);
4203 neighbour_node_number = i0;
4205 neighbour_node_number, value_id);
4209 neighbour_node_number = neigh_n_p - 1 + neigh_n_p * i0;
4211 neighbour_node_number, value_id);
4215 neighbour_node_number = neigh_n_p * i0;
4217 neighbour_node_number, value_id);
4222 OOMPH_CURRENT_FUNCTION,
4223 OOMPH_EXCEPTION_LOCATION);
4227 master_node_number.push_back(neighbour_node_number);
4228 master_node_pt.push_back(neighbour_node_pt);
4232 unsigned local_node_number = 0;
4233 Node* local_node_pt = 0;
4236 for (
unsigned i0 = 0; i0 < my_n_p; i0++)
4247 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
4248 s_fraction[1] = 1.0;
4249 local_node_number = i0 + my_n_p * (my_n_p - 1);
4251 this->interpolating_node_pt(local_node_number, value_id);
4256 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
4257 s_fraction[1] = 0.0;
4258 local_node_number = i0;
4260 this->interpolating_node_pt(local_node_number, value_id);
4264 s_fraction[0] = 1.0;
4266 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
4267 local_node_number = my_n_p - 1 + my_n_p * i0;
4269 this->interpolating_node_pt(local_node_number, value_id);
4273 s_fraction[0] = 0.0;
4275 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
4276 local_node_number = my_n_p * i0;
4278 this->interpolating_node_pt(local_node_number, value_id);
4283 OOMPH_CURRENT_FUNCTION,
4284 OOMPH_EXCEPTION_LOCATION);
4288 dependent_node_number.push_back(local_node_number);
4289 dependent_node_pt.push_back(local_node_pt);
4292 dependent_node_s_fraction.push_back(s_fraction);
4296 const unsigned n_dependent_nodes = dependent_node_pt.size();
4297 const unsigned n_master_nodes = master_node_pt.size();
4298 const unsigned dependent_element_nnode_1d = my_n_p;
4299 const unsigned master_element_nnode_1d = neigh_n_p;
4309 dependent_element_nnode_1d, dependent_nodal_position, dependent_weight);
4313 master_element_nnode_1d, master_nodal_position, master_weight);
4326 const unsigned n_mortar_vertices = 2;
4329 vertex_pos[1] = this->ninterpolating_node_1d(value_id) - 1;
4331 for (
unsigned i = 0;
i < my_n_p - n_mortar_vertices;
i++)
4333 non_vertex_pos[
i] =
i + 1;
4337 for (
unsigned v = 0; v < n_mortar_vertices; v++)
4340 bool node_is_shared =
false;
4341 for (
unsigned i = 0;
i < master_node_pt.size();
i++)
4343 if (dependent_node_pt[vertex_pos[v]] == master_node_pt[
i])
4345 node_is_shared =
true;
4352 if (!node_is_shared)
4359 for (
unsigned i = 0;
i < 2;
i++)
4363 dependent_node_s_fraction[vertex_pos[v]][translate_s[
i]] *
4364 (s_hi_neigh[
i] - s_lo_neigh[
i]);
4369 s_in_neigh, master_shapes, value_id);
4373 dependent_node_pt[vertex_pos[v]]->set_nonhanging();
4377 for (
unsigned m = 0; m < n_master_nodes; m++)
4380 m, master_node_pt[m], master_shapes[master_node_number[m]]);
4384 dependent_node_pt[vertex_pos[v]]->set_hanging_pt(explicit_hang_pt,
4391 if (n_dependent_nodes - n_mortar_vertices > 0)
4398 for (
unsigned i = 0;
i < shared_node_M.size();
i++)
4400 shared_node_M[
i].resize(n_dependent_nodes - n_mortar_vertices);
4404 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
4408 pow(-1.0,
int((dependent_element_nnode_1d - 1) -
i - 1)) *
4410 dependent_nodal_position[non_vertex_pos[
i]]);
4412 diag_M[
i] = psi[
i] * dependent_weight[non_vertex_pos[
i]];
4416 for (
unsigned v = 0; v < shared_node_M.size(); v++)
4418 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
4421 if (std::fabs(dependent_nodal_position[non_vertex_pos[
i]] -
4422 dependent_nodal_position[vertex_pos[v]]) >= 1.0e-8)
4426 pow(-1.0,
int((dependent_element_nnode_1d - 1) -
i - 1)) *
4428 dependent_nodal_position[vertex_pos[v]]) /
4429 (dependent_nodal_position[non_vertex_pos[
i]] -
4430 dependent_nodal_position[vertex_pos[v]]);
4434 dependent_element_nnode_1d - 1,
4435 dependent_nodal_position[vertex_pos[v]])) < 1.0e-8)
4439 pow(-1.0,
int((dependent_element_nnode_1d - 1) -
i - 1)) *
4441 dependent_element_nnode_1d - 1,
4442 dependent_nodal_position[non_vertex_pos[
i]]);
4448 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4449 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4450 OOMPH_EXCEPTION_LOCATION);
4453 shared_node_M[v][
i] = psi[
i] * dependent_weight[vertex_pos[v]];
4463 for (
unsigned i = 0;
i < P.size();
i++)
4465 P[
i].resize(n_master_nodes, 0.0);
4480 if (dependent_element_nnode_1d >= master_element_nnode_1d)
4483 quadrature_knot = &dependent_nodal_position;
4484 quadrature_weight = &dependent_weight;
4489 quadrature_knot = &master_nodal_position;
4490 quadrature_weight = &master_weight;
4494 for (
unsigned q = 0; q < (*quadrature_weight).size(); q++)
4500 s_on_mortar[0] = (*quadrature_knot)[q];
4503 for (
unsigned k = 0; k < n_dependent_nodes - n_mortar_vertices; k++)
4506 if (std::fabs(dependent_nodal_position[non_vertex_pos[k]] -
4507 s_on_mortar[0]) >= 1.0e-08)
4511 pow(-1.0,
int((dependent_element_nnode_1d - 1) - k - 1)) *
4514 (dependent_nodal_position[non_vertex_pos[k]] - s_on_mortar[0]);
4518 dependent_element_nnode_1d - 1, s_on_mortar[0])) <
4523 pow(-1.0,
int((dependent_element_nnode_1d - 1) - k - 1)) *
4531 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4532 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4533 OOMPH_EXCEPTION_LOCATION);
4539 for (
unsigned i = 0;
i < 2;
i++)
4541 s_fraction[
i] = (
i == active_coord_index) ?
4542 0.5 * (s_on_mortar[0] + 1.0) :
4543 dependent_node_s_fraction[vertex_pos[0]][
i];
4548 for (
unsigned i = 0;
i < 2;
i++)
4550 s_in_neigh[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]] *
4551 (s_hi_neigh[
i] - s_lo_neigh[
i]);
4556 s_in_neigh, master_shapes, value_id);
4559 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
4561 for (
unsigned j = 0; j < n_master_nodes; j++)
4563 P[
i][j] += master_shapes[master_node_number[j]] * psi[
i] *
4564 (*quadrature_weight)[q];
4575 for (
unsigned v = 0; v < n_mortar_vertices; v++)
4579 for (
unsigned i = 0;
i < 2;
i++)
4582 (
i == active_coord_index) ?
4583 0.5 * (dependent_nodal_position[vertex_pos[v]] + 1.0) :
4584 dependent_node_s_fraction[vertex_pos[0]][
i];
4589 for (
unsigned i = 0;
i < 2;
i++)
4591 s_in_neigh[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]] *
4592 (s_hi_neigh[
i] - s_lo_neigh[
i]);
4597 s_in_neigh, master_shapes, value_id);
4599 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
4601 for (
unsigned k = 0; k < n_master_nodes; k++)
4604 master_shapes[master_node_number[k]] * shared_node_M[v][
i];
4611 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
4613 for (
unsigned j = 0; j < n_master_nodes; j++)
4615 P[
i][j] /= diag_M[
i];
4622 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
4624 hang_info_pt[
i] =
new HangInfo(n_master_nodes);
4629 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
4631 for (
unsigned j = 0; j < n_master_nodes; j++)
4633 hang_info_pt[
i]->set_master_node_pt(j, master_node_pt[j], P[
i][j]);
4639 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
4645 bool node_is_really_shared =
false;
4646 for (
unsigned m = 0; m < hang_info_pt[
i]->nmaster(); m++)
4650 if (hang_info_pt[
i]->master_node_pt(m) ==
4651 dependent_node_pt[non_vertex_pos[
i]])
4653 node_is_really_shared =
true;
4657 if (std::fabs(hang_info_pt[
i]->master_weight(m) - 1.0) > 1.0e-06)
4660 "Something fishy here -- with shared node at a mortar vertex",
4661 "PRefineableQElemen<2,INITIAL_NNODE_1D>t::quad_hang_helper()",
4662 OOMPH_EXCEPTION_LOCATION);
4669 if (!node_is_really_shared)
4671 dependent_node_pt[non_vertex_pos[
i]]->set_hanging_pt(
4672 hang_info_pt[
i], -1);
4680 for (
unsigned i = 0;
i < n_dependent_nodes;
i++)
4683 if (dependent_node_pt[
i]->is_hanging())
4686 double weight_sum = 0.0;
4687 for (
unsigned m = 0;
4688 m < dependent_node_pt[
i]->hanging_pt()->nmaster();
4691 weight_sum += dependent_node_pt[
i]->hanging_pt()->master_weight(m);
4695 if (std::fabs(weight_sum - 1.0) > 1.0e-08)
4697 oomph_info <<
"Sum of master weights: " << weight_sum << std::endl;
4699 "Weights in hanging scheme do not sum to 1",
4700 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4701 OOMPH_EXCEPTION_LOCATION);
4708 bool is_master =
false;
4709 for (
unsigned n = 0; n < n_master_nodes; n++)
4711 if (dependent_node_pt[
i] == master_node_pt[n])
4722 std::ostringstream error_string;
4723 error_string <<
"This node in the dependent element is neither"
4725 <<
"hanging or shared with a master element."
4730 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4731 OOMPH_EXCEPTION_LOCATION);
4744 for (
unsigned i = 0;
i < n_dependent_nodes;
i++)
4747 if (dependent_node_pt[
i]->is_hanging())
4754 this->local_coordinate_of_node(dependent_node_number[
i], s_local);
4759 this->interpolated_x(s_local, x_in_neighb);
4763 dependent_node_pt[
i]->x(0) = x_in_neighb[0];
4764 dependent_node_pt[
i]->x(1) = x_in_neighb[1];
4776 template<
unsigned INITIAL_NNODE_1D>
4781 unsigned Nnode_1d = this->nnode_1d();
4782 unsigned n0 = n % Nnode_1d;
4783 unsigned n1 = unsigned(
double(n) /
double(Nnode_1d)) % Nnode_1d;
4784 unsigned n2 = unsigned(
double(n) /
double(Nnode_1d * Nnode_1d));
4825 std::ostringstream error_message;
4826 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
4827 error_message <<
"\n shape functions.\n";
4829 OOMPH_CURRENT_FUNCTION,
4830 OOMPH_EXCEPTION_LOCATION);
4836 template<
unsigned INITIAL_NNODE_1D>
4840 this->local_coordinate_of_node(n, s_fraction);
4841 s_fraction[0] = 0.5 * (s_fraction[0] + 1.0);
4842 s_fraction[1] = 0.5 * (s_fraction[1] + 1.0);
4843 s_fraction[2] = 0.5 * (s_fraction[2] + 1.0);
4846 template<
unsigned INITIAL_NNODE_1D>
4848 const unsigned& n1d,
const unsigned&
i)
4850 switch (this->nnode_1d())
4877 std::ostringstream error_message;
4878 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
4879 error_message <<
"\n shape functions.\n";
4881 OOMPH_CURRENT_FUNCTION,
4882 OOMPH_EXCEPTION_LOCATION);
4890 template<
unsigned INITIAL_NNODE_1D>
4894 unsigned Nnode_1d = this->nnode_1d();
4901 for (
unsigned i = 0;
i < 3;
i++)
4906 bool is_found =
false;
4909 if (std::fabs(
s[
i] + 1.0) < tol)
4915 else if (std::fabs(
s[
i] - 1.0) < tol)
4917 index[
i] = Nnode_1d - 1;
4927 for (
unsigned n = 1; n < Nnode_1d - 1; n++)
4929 if (std::fabs(z[n] -
s[
i]) < tol)
4945 return this->node_pt(index[0] + Nnode_1d * index[1] +
4946 Nnode_1d * Nnode_1d * index[2]);
4955 template<
unsigned INITIAL_NNODE_1D>
4959 using namespace OcTreeNames;
4965 if (s_fraction[0] == 0.0)
4968 if (s_fraction[1] == 0.0)
4970 edges.push_back(
LD);
4972 if (s_fraction[2] == 0.0)
4974 edges.push_back(
LB);
4976 if (s_fraction[1] == 1.0)
4978 edges.push_back(
LU);
4980 if (s_fraction[2] == 1.0)
4982 edges.push_back(
LF);
4986 if (s_fraction[0] == 1.0)
4989 if (s_fraction[1] == 0.0)
4991 edges.push_back(
RD);
4993 if (s_fraction[2] == 0.0)
4995 edges.push_back(
RB);
4997 if (s_fraction[1] == 1.0)
4999 edges.push_back(
RU);
5001 if (s_fraction[2] == 1.0)
5003 edges.push_back(
RF);
5007 if (s_fraction[1] == 0.0)
5010 if (s_fraction[2] == 0.0)
5012 edges.push_back(
DB);
5014 if (s_fraction[2] == 1.0)
5016 edges.push_back(
DF);
5020 if (s_fraction[1] == 1.0)
5023 if (s_fraction[2] == 0.0)
5025 edges.push_back(
UB);
5027 if (s_fraction[2] == 1.0)
5029 edges.push_back(
UF);
5033 if (s_fraction[2] == 0.0)
5038 if (s_fraction[2] == 1.0)
5044 unsigned n_face = faces.size();
5047 unsigned n_edge = edges.size();
5054 int neigh_face, diff_level;
5055 bool in_neighbouring_tree;
5059 for (
unsigned j = 0; j < n_face; j++)
5069 in_neighbouring_tree);
5078 bool a_vertex_node_is_built =
false;
5081 if (neigh_obj_pt == 0)
5084 OOMPH_CURRENT_FUNCTION,
5085 OOMPH_EXCEPTION_LOCATION);
5087 for (
unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node(); vnode++)
5089 if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
5090 a_vertex_node_is_built =
true;
5093 if (a_vertex_node_is_built)
5103 for (
unsigned i = 0;
i < 3;
i++)
5105 s[
i] = s_lo_neigh[
i] +
5106 s_fraction[translate_s[
i]] * (s_hi_neigh[
i] - s_lo_neigh[
i]);
5110 Node* neighbour_node_pt =
5114 if (neighbour_node_pt != 0)
5118 if (in_neighbouring_tree)
5122 octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
5125 return neighbour_node_pt;
5134 for (
unsigned j = 0; j < n_edge; j++)
5148 unsigned i_root_edge_neighbour = 0;
5151 unsigned nroot_edge_neighbour = 0;
5155 bool keep_searching =
true;
5156 while (keep_searching)
5161 i_root_edge_neighbour,
5162 nroot_edge_neighbour,
5176 bool a_vertex_node_is_built =
false;
5179 if (neigh_obj_pt == 0)
5182 OOMPH_CURRENT_FUNCTION,
5183 OOMPH_EXCEPTION_LOCATION);
5185 for (
unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node();
5188 if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
5189 a_vertex_node_is_built =
true;
5192 if (a_vertex_node_is_built)
5202 for (
unsigned i = 0;
i < 3;
i++)
5204 s[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]] *
5205 (s_hi_neigh[
i] - s_lo_neigh[
i]);
5209 Node* neighbour_node_pt =
5213 if (neighbour_node_pt != 0)
5220 unsigned n_faces_attached_to_edge = faces_attached_to_edge.size();
5223 for (
unsigned i_face = 0; i_face < n_faces_attached_to_edge;
5227 is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
5228 faces_attached_to_edge[i_face]);
5239 return neighbour_node_pt;
5246 i_root_edge_neighbour++;
5249 if (i_root_edge_neighbour >= nroot_edge_neighbour)
5251 keep_searching =
false;
5269 template<
unsigned INITIAL_NNODE_1D>
5274 using namespace OcTreeNames;
5280 if (s_fraction[0] == 0.0)
5283 if (s_fraction[1] == 0.0)
5285 edges.push_back(
LD);
5287 if (s_fraction[2] == 0.0)
5289 edges.push_back(
LB);
5291 if (s_fraction[1] == 1.0)
5293 edges.push_back(
LU);
5295 if (s_fraction[2] == 1.0)
5297 edges.push_back(
LF);
5301 if (s_fraction[0] == 1.0)
5304 if (s_fraction[1] == 0.0)
5306 edges.push_back(
RD);
5308 if (s_fraction[2] == 0.0)
5310 edges.push_back(
RB);
5312 if (s_fraction[1] == 1.0)
5314 edges.push_back(
RU);
5316 if (s_fraction[2] == 1.0)
5318 edges.push_back(
RF);
5322 if (s_fraction[1] == 0.0)
5325 if (s_fraction[2] == 0.0)
5327 edges.push_back(
DB);
5329 if (s_fraction[2] == 1.0)
5331 edges.push_back(
DF);
5335 if (s_fraction[1] == 1.0)
5338 if (s_fraction[2] == 0.0)
5340 edges.push_back(
UB);
5342 if (s_fraction[2] == 1.0)
5344 edges.push_back(
UF);
5348 if (s_fraction[2] == 0.0)
5353 if (s_fraction[2] == 1.0)
5359 unsigned n_face = faces.size();
5362 unsigned n_edge = edges.size();
5369 int neigh_face, diff_level;
5370 bool in_neighbouring_tree;
5374 for (
unsigned j = 0; j < n_face; j++)
5384 in_neighbouring_tree);
5401 for (
unsigned i = 0;
i < 3;
i++)
5403 s[
i] = s_lo_neigh[
i] +
5404 s_fraction[translate_s[
i]] * (s_hi_neigh[
i] - s_lo_neigh[
i]);
5408 if (neigh_pt->
nsons() != 0)
5426 s_in_son[0] = 1.0 + 2.0 *
s[0];
5427 s_in_son[1] = 1.0 + 2.0 *
s[1];
5428 s_in_son[2] = 1.0 + 2.0 *
s[2];
5435 s_in_son[0] = 1.0 + 2.0 *
s[0];
5436 s_in_son[1] = 1.0 + 2.0 *
s[1];
5437 s_in_son[2] = -1.0 + 2.0 *
s[2];
5448 s_in_son[0] = 1.0 + 2.0 *
s[0];
5449 s_in_son[1] = -1.0 + 2.0 *
s[1];
5450 s_in_son[2] = 1.0 + 2.0 *
s[2];
5457 s_in_son[0] = 1.0 + 2.0 *
s[0];
5458 s_in_son[1] = -1.0 + 2.0 *
s[1];
5459 s_in_son[2] = -1.0 + 2.0 *
s[2];
5474 s_in_son[0] = -1.0 + 2.0 *
s[0];
5475 s_in_son[1] = 1.0 + 2.0 *
s[1];
5476 s_in_son[2] = 1.0 + 2.0 *
s[2];
5483 s_in_son[0] = -1.0 + 2.0 *
s[0];
5484 s_in_son[1] = 1.0 + 2.0 *
s[1];
5485 s_in_son[2] = -1.0 + 2.0 *
s[2];
5496 s_in_son[0] = -1.0 + 2.0 *
s[0];
5497 s_in_son[1] = -1.0 + 2.0 *
s[1];
5498 s_in_son[2] = 1.0 + 2.0 *
s[2];
5505 s_in_son[0] = -1.0 + 2.0 *
s[0];
5506 s_in_son[1] = -1.0 + 2.0 *
s[1];
5507 s_in_son[2] = -1.0 + 2.0 *
s[2];
5513 Node* neighbour_son_node_pt =
5518 if (neighbour_son_node_pt != 0)
5522 if (in_neighbouring_tree)
5526 octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
5530 return neighbour_son_node_pt;
5540 for (
unsigned j = 0; j < n_edge; j++)
5554 unsigned i_root_edge_neighbour = 0;
5557 unsigned nroot_edge_neighbour = 0;
5561 bool keep_searching =
true;
5562 while (keep_searching)
5567 i_root_edge_neighbour,
5568 nroot_edge_neighbour,
5590 for (
unsigned i = 0;
i < 3;
i++)
5592 s[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]] *
5593 (s_hi_neigh[
i] - s_lo_neigh[
i]);
5597 if (neigh_pt->
nsons() != 0)
5615 s_in_son[0] = 1.0 + 2.0 *
s[0];
5616 s_in_son[1] = 1.0 + 2.0 *
s[1];
5617 s_in_son[2] = 1.0 + 2.0 *
s[2];
5624 s_in_son[0] = 1.0 + 2.0 *
s[0];
5625 s_in_son[1] = 1.0 + 2.0 *
s[1];
5626 s_in_son[2] = -1.0 + 2.0 *
s[2];
5637 s_in_son[0] = 1.0 + 2.0 *
s[0];
5638 s_in_son[1] = -1.0 + 2.0 *
s[1];
5639 s_in_son[2] = 1.0 + 2.0 *
s[2];
5646 s_in_son[0] = 1.0 + 2.0 *
s[0];
5647 s_in_son[1] = -1.0 + 2.0 *
s[1];
5648 s_in_son[2] = -1.0 + 2.0 *
s[2];
5663 s_in_son[0] = -1.0 + 2.0 *
s[0];
5664 s_in_son[1] = 1.0 + 2.0 *
s[1];
5665 s_in_son[2] = 1.0 + 2.0 *
s[2];
5672 s_in_son[0] = -1.0 + 2.0 *
s[0];
5673 s_in_son[1] = 1.0 + 2.0 *
s[1];
5674 s_in_son[2] = -1.0 + 2.0 *
s[2];
5685 s_in_son[0] = -1.0 + 2.0 *
s[0];
5686 s_in_son[1] = -1.0 + 2.0 *
s[1];
5687 s_in_son[2] = 1.0 + 2.0 *
s[2];
5694 s_in_son[0] = -1.0 + 2.0 *
s[0];
5695 s_in_son[1] = -1.0 + 2.0 *
s[1];
5696 s_in_son[2] = -1.0 + 2.0 *
s[2];
5702 Node* neighbour_son_node_pt =
5708 if (neighbour_son_node_pt != 0)
5715 unsigned n_faces_attached_to_edge =
5716 faces_attached_to_edge.size();
5719 for (
unsigned i_face = 0; i_face < n_faces_attached_to_edge;
5723 is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
5724 faces_attached_to_edge[i_face]);
5737 return neighbour_son_node_pt;
5745 i_root_edge_neighbour++;
5748 if (i_root_edge_neighbour >= nroot_edge_neighbour)
5750 keep_searching =
false;
5767 template<
unsigned INITIAL_NNODE_1D>
5769 Tree*
const& adopted_father_pt,
const unsigned& initial_p_order)
5775 if (adopted_father_pt != 0)
5778 father_pt =
dynamic_cast<OcTree*
>(adopted_father_pt);
5781 else if (Tree_pt != 0)
5784 father_pt =
dynamic_cast<OcTree*
>(octree_pt()->father_pt());
5796 if (father_pt != 0 || initial_p_order != 0)
5803 if (father_el_pt != 0)
5805 unsigned father_p_order = father_el_pt->
p_order();
5807 P_order = father_p_order;
5812 P_order = initial_p_order;
5817 unsigned new_n_node = P_order * P_order * P_order;
5820 this->set_n_node(new_n_node);
5823 delete this->integral_pt();
5845 std::ostringstream error_message;
5846 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
5847 error_message <<
"\n integration scheme.\n";
5849 OOMPH_CURRENT_FUNCTION,
5850 OOMPH_EXCEPTION_LOCATION);
5859 template<
unsigned INITIAL_NNODE_1D>
5863 using namespace OcTreeNames;
5866 unsigned n_p = nnode_1d();
5869 if (Father_bound[n_p].nrow() == 0)
5871 setup_father_bounds();
5875 OcTree* father_pt =
dynamic_cast<OcTree*
>(octree_pt()->father_pt());
5878 int son_type = Tree_pt->
son_type();
5887 "Something fishy here: I have no father and yet \n";
5888 error_message +=
"I have no nodes. Who has created me then?!\n";
5892 "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
5893 OOMPH_EXCEPTION_LOCATION);
5907 unsigned ntstorage = time_stepper_pt->ntstorage();
6013 if (father_el_pt->
node_pt(0) == 0)
6016 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
6017 "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
6018 OOMPH_EXCEPTION_LOCATION);
6026 for (
unsigned i0 = 0; i0 < n_p; i0++)
6029 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
6031 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
6033 for (
unsigned i1 = 0; i1 < n_p; i1++)
6036 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
6038 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
6040 for (
unsigned i2 = 0; i2 < n_p; i2++)
6044 s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
6046 s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
6049 jnod = i0 + n_p * i1 + n_p * n_p * i2;
6053 Node* created_node_pt =
6058 if (created_node_pt != 0)
6061 this->node_pt(jnod) = created_node_pt;
6069 for (
unsigned t = 0;
t < ntstorage;
t++)
6078 unsigned n_val_at_node = created_node_pt->
nvalue();
6079 unsigned n_val_from_function = prev_values.size();
6081 unsigned n_var = n_val_at_node < n_val_from_function ?
6083 n_val_from_function;
6085 for (
unsigned k = 0; k < n_var; k++)
6087 created_node_pt->
set_value(
t, k, prev_values[k]);
6103 template<
unsigned INITIAL_NNODE_1D>
6115 if (clone_el_pt == 0)
6118 "Cloned copy must be a PRefineableQElement<3,INITIAL_NNODE_1D>!",
6119 OOMPH_CURRENT_FUNCTION,
6120 OOMPH_EXCEPTION_LOCATION);
6127 bool clone_is_ok =
true;
6130 clone_is_ok = clone_is_ok && (clone_el_pt->
p_order() == this->p_order());
6134 std::ostringstream err_stream;
6135 err_stream <<
"Clone element has a different p-order from me,"
6137 <<
"but it is supposed to be a copy!" << std::endl;
6140 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6144 clone_is_ok = clone_is_ok && (clone_el_pt->
nnode() == this->nnode());
6148 std::ostringstream err_stream;
6149 err_stream <<
"Clone element has a different number of nodes from me,"
6151 <<
"but it is supposed to be a copy!" << std::endl;
6154 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6158 for (
unsigned n = 0; n < this->nnode(); n++)
6161 clone_is_ok && (clone_el_pt->
node_pt(n) == this->node_pt(n));
6166 std::ostringstream err_stream;
6167 err_stream <<
"The nodes belonging to the clone element are different"
6169 <<
"from mine, but it is supposed to be a copy!"
6173 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6185 throw OomphLibError(
"Can't handle generalised nodal positions (yet).",
6186 OOMPH_CURRENT_FUNCTION,
6187 OOMPH_EXCEPTION_LOCATION);
6192 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
6195 unsigned ntstorage = time_stepper_pt->
ntstorage();
6198 this->P_order += inc;
6201 delete this->integral_pt();
6202 switch (this->P_order)
6223 std::ostringstream error_message;
6224 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
6225 error_message <<
"\n integration scheme.\n";
6227 OOMPH_CURRENT_FUNCTION,
6228 OOMPH_EXCEPTION_LOCATION);
6232 this->set_n_node(this->P_order * this->P_order * this->P_order);
6255 for (
unsigned i0 = 0; i0 < this->P_order; i0++)
6258 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
6260 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
6262 for (
unsigned i1 = 0; i1 < this->P_order; i1++)
6265 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
6267 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
6269 for (
unsigned i2 = 0; i2 < this->P_order; i2++)
6272 s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
6274 s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
6277 jnod = i0 + this->P_order * i1 + this->P_order * this->P_order * i2;
6281 bool node_done =
false;
6290 if (created_node_pt != 0)
6293 this->node_pt(jnod) = created_node_pt;
6301 for (
unsigned t = 0;
t < ntstorage;
t++)
6310 unsigned n_val_at_node = created_node_pt->
nvalue();
6311 unsigned n_val_from_function = prev_values.size();
6313 unsigned n_var = n_val_at_node < n_val_from_function ?
6315 n_val_from_function;
6317 for (
unsigned k = 0; k < n_var; k++)
6319 created_node_pt->
set_value(
t, k, prev_values[k]);
6335 bool is_periodic =
false;
6337 this->node_created_by_neighbour(s_fraction, is_periodic);
6340 if (created_node_pt != 0)
6349 Node* neighbour_node_pt = created_node_pt;
6360 else if (i0 == this->P_order - 1)
6385 else if (i1 == this->P_order - 1)
6440 else if (i2 == this->P_order - 1)
6479 std::set<unsigned> boundaries;
6491 if (boundaries.size() > 2)
6494 "boundaries.size()>2 seems a bit strange..\n",
6495 OOMPH_CURRENT_FUNCTION,
6496 OOMPH_EXCEPTION_LOCATION);
6500 if (boundaries.size() == 0)
6502 std::ostringstream error_stream;
6503 error_stream <<
"Periodic node is not on a boundary...\n"
6504 <<
"Coordinates: " << created_node_pt->
x(0)
6505 <<
" " << created_node_pt->
x(1) <<
"\n";
6507 OOMPH_CURRENT_FUNCTION,
6508 OOMPH_EXCEPTION_LOCATION);
6514 this->construct_boundary_node(jnod, time_stepper_pt);
6519 for (
unsigned t = 0;
t < ntstorage;
t++)
6527 clone_el_pt->
get_x(
t,
s, x_prev);
6529 for (
unsigned i = 0;
i < n_dim;
i++)
6531 created_node_pt->
x(
t,
i) = x_prev[
i];
6540 for (std::set<unsigned>::iterator it = boundaries.begin();
6541 it != boundaries.end();
6555 *it, my_bound,
s, zeta);
6569 this->node_pt(jnod) = created_node_pt;
6583 bool is_periodic =
false;
6585 this->node_created_by_son_of_neighbour(s_fraction, is_periodic);
6588 if (created_node_pt != 0)
6597 Node* neighbour_node_pt = created_node_pt;
6608 else if (i0 == this->P_order - 1)
6633 else if (i1 == this->P_order - 1)
6688 else if (i2 == this->P_order - 1)
6727 std::set<unsigned> boundaries;
6739 if (boundaries.size() > 2)
6742 "boundaries.size()>2 seems a bit strange..\n",
6743 OOMPH_CURRENT_FUNCTION,
6744 OOMPH_EXCEPTION_LOCATION);
6748 if (boundaries.size() == 0)
6750 std::ostringstream error_stream;
6751 error_stream <<
"Periodic node is not on a boundary...\n"
6752 <<
"Coordinates: " << created_node_pt->
x(0)
6753 <<
" " << created_node_pt->
x(1) <<
"\n";
6755 OOMPH_CURRENT_FUNCTION,
6756 OOMPH_EXCEPTION_LOCATION);
6762 this->construct_boundary_node(jnod, time_stepper_pt);
6767 for (
unsigned t = 0;
t < ntstorage;
t++)
6775 clone_el_pt->
get_x(
t,
s, x_prev);
6777 for (
unsigned i = 0;
i < n_dim;
i++)
6779 created_node_pt->
x(
t,
i) = x_prev[
i];
6788 for (std::set<unsigned>::iterator it = boundaries.begin();
6789 it != boundaries.end();
6803 *it, my_bound,
s, zeta);
6817 this->node_pt(jnod) = created_node_pt;
6840 else if (i0 == this->P_order - 1)
6864 else if (i1 == this->P_order - 1)
6918 else if (i2 == this->P_order - 1)
6957 std::set<unsigned> boundaries;
6967 if (boundaries.size() > 2)
6970 "boundaries.size()>2 seems a bit strange..\n",
6971 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6972 OOMPH_EXCEPTION_LOCATION);
6978 if (boundaries.size() > 0)
6982 this->construct_boundary_node(jnod, time_stepper_pt);
6989 Vector<int> bound_cons(this->ncont_interpolated_values());
6990 clone_el_pt->
get_bcs(my_bound, bound_cons);
6993 unsigned n_value = created_node_pt->
nvalue();
6994 for (
unsigned k = 0; k < n_value; k++)
6998 created_node_pt->
pin(k);
7005 dynamic_cast<SolidNode*
>(created_node_pt);
7006 if (solid_node_pt != 0)
7009 unsigned n_dim = created_node_pt->
ndim();
7014 if (clone_solid_el_pt == 0)
7017 "We have a SolidNode outside a refineable SolidElement\n";
7019 "during mesh refinement -- this doesn't make sense";
7023 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
7024 OOMPH_EXCEPTION_LOCATION);
7027 clone_solid_el_pt->
get_solid_bcs(my_bound, solid_bound_cons);
7030 for (
unsigned k = 0; k < n_dim; k++)
7032 if (solid_bound_cons[k])
7045 for (std::set<unsigned>::iterator it = boundaries.begin();
7046 it != boundaries.end();
7060 *it, my_bound,
s, zeta);
7072 created_node_pt = this->construct_node(jnod, time_stepper_pt);
7088 for (
unsigned t = 0;
t < ntstorage;
t++)
7096 clone_el_pt->
get_x(
t,
s, x_prev);
7099 for (
unsigned i = 0;
i < 3;
i++)
7101 created_node_pt->
x(
t,
i) = x_prev[
i];
7106 for (
unsigned t = 0;
t < ntstorage;
t++)
7114 unsigned n_value = created_node_pt->
nvalue();
7115 for (
unsigned k = 0; k < n_value; k++)
7117 created_node_pt->
set_value(
t, k, prev_values[k]);
7134 "Have not implemented p-refinement for";
7135 error_message +=
"Algebraic p-refineable elements yet\n";
7139 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
7140 OOMPH_EXCEPTION_LOCATION);
7164 this->node_pt(jnod),
s, clone_el_pt);
7178 for (
unsigned i0 = 0; i0 < this->P_order; i0++)
7181 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
7183 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
7185 for (
unsigned i1 = 0; i1 < this->P_order; i1++)
7188 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
7190 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
7192 for (
unsigned i2 = 0; i2 < this->P_order; i2++)
7195 s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
7197 s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
7200 jnod = i0 + this->P_order * i1 + this->P_order * this->P_order * i2;
7203 for (
unsigned t = 0;
t < ntstorage;
t++)
7211 this->get_x(
t,
s, x_prev);
7214 for (
unsigned i = 0;
i < 3;
i++)
7216 this->node_pt(jnod)->x(
t,
i) = x_prev[
i];
7232 if (clone_m_el_pt != 0)
7246 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
7248 "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
7249 error_message +=
"then I should be too....\n";
7253 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
7254 OOMPH_EXCEPTION_LOCATION);
7272 this->further_build();
7278 template<
unsigned INITIAL_NNODE_1D>
7292 for (
unsigned i = 0;
i < P_order;
i++)
7294 for (
unsigned j = 0; j < P_order; j++)
7296 for (
unsigned k = 0; k < P_order; k++)
7298 psi(P_order * P_order *
i + P_order * j + k) =
7299 psi3[
i] * psi2[j] * psi1[k];
7313 for (
unsigned i = 0;
i < P_order;
i++)
7315 for (
unsigned j = 0; j < P_order; j++)
7317 for (
unsigned k = 0; k < P_order; k++)
7319 psi(P_order * P_order *
i + P_order * j + k) =
7320 psi3[
i] * psi2[j] * psi1[k];
7334 for (
unsigned i = 0;
i < P_order;
i++)
7336 for (
unsigned j = 0; j < P_order; j++)
7338 for (
unsigned k = 0; k < P_order; k++)
7340 psi(P_order * P_order *
i + P_order * j + k) =
7341 psi3[
i] * psi2[j] * psi1[k];
7355 for (
unsigned i = 0;
i < P_order;
i++)
7357 for (
unsigned j = 0; j < P_order; j++)
7359 for (
unsigned k = 0; k < P_order; k++)
7361 psi(P_order * P_order *
i + P_order * j + k) =
7362 psi3[
i] * psi2[j] * psi1[k];
7376 for (
unsigned i = 0;
i < P_order;
i++)
7378 for (
unsigned j = 0; j < P_order; j++)
7380 for (
unsigned k = 0; k < P_order; k++)
7382 psi(P_order * P_order *
i + P_order * j + k) =
7383 psi3[
i] * psi2[j] * psi1[k];
7397 for (
unsigned i = 0;
i < P_order;
i++)
7399 for (
unsigned j = 0; j < P_order; j++)
7401 for (
unsigned k = 0; k < P_order; k++)
7403 psi(P_order * P_order *
i + P_order * j + k) =
7404 psi3[
i] * psi2[j] * psi1[k];
7411 std::ostringstream error_message;
7412 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
7413 error_message <<
"\n polynomial order for shape functions.\n";
7415 OOMPH_CURRENT_FUNCTION,
7416 OOMPH_EXCEPTION_LOCATION);
7423 template<
unsigned INITIAL_NNODE_1D>
7440 for (
unsigned i = 0;
i < P_order;
i++)
7442 for (
unsigned j = 0; j < P_order; j++)
7444 for (
unsigned k = 0; k < P_order; k++)
7447 dpsids(index, 0) = psi3[
i] * psi2[j] * dpsi1ds[k];
7448 dpsids(index, 1) = psi3[
i] * dpsi2ds[j] * psi1[k];
7449 dpsids(index, 2) = dpsi3ds[
i] * psi2[j] * psi1[k];
7450 psi[index] = psi3[
i] * psi2[j] * psi1[k];
7469 for (
unsigned i = 0;
i < P_order;
i++)
7471 for (
unsigned j = 0; j < P_order; j++)
7473 for (
unsigned k = 0; k < P_order; k++)
7476 dpsids(index, 0) = psi3[
i] * psi2[j] * dpsi1ds[k];
7477 dpsids(index, 1) = psi3[
i] * dpsi2ds[j] * psi1[k];
7478 dpsids(index, 2) = dpsi3ds[
i] * psi2[j] * psi1[k];
7479 psi[index] = psi3[
i] * psi2[j] * psi1[k];
7498 for (
unsigned i = 0;
i < P_order;
i++)
7500 for (
unsigned j = 0; j < P_order; j++)
7502 for (
unsigned k = 0; k < P_order; k++)
7505 dpsids(index, 0) = psi3[
i] * psi2[j] * dpsi1ds[k];
7506 dpsids(index, 1) = psi3[
i] * dpsi2ds[j] * psi1[k];
7507 dpsids(index, 2) = dpsi3ds[
i] * psi2[j] * psi1[k];
7508 psi[index] = psi3[
i] * psi2[j] * psi1[k];
7527 for (
unsigned i = 0;
i < P_order;
i++)
7529 for (
unsigned j = 0; j < P_order; j++)
7531 for (
unsigned k = 0; k < P_order; k++)
7534 dpsids(index, 0) = psi3[
i] * psi2[j] * dpsi1ds[k];
7535 dpsids(index, 1) = psi3[
i] * dpsi2ds[j] * psi1[k];
7536 dpsids(index, 2) = dpsi3ds[
i] * psi2[j] * psi1[k];
7537 psi[index] = psi3[
i] * psi2[j] * psi1[k];
7556 for (
unsigned i = 0;
i < P_order;
i++)
7558 for (
unsigned j = 0; j < P_order; j++)
7560 for (
unsigned k = 0; k < P_order; k++)
7563 dpsids(index, 0) = psi3[
i] * psi2[j] * dpsi1ds[k];
7564 dpsids(index, 1) = psi3[
i] * dpsi2ds[j] * psi1[k];
7565 dpsids(index, 2) = dpsi3ds[
i] * psi2[j] * psi1[k];
7566 psi[index] = psi3[
i] * psi2[j] * psi1[k];
7585 for (
unsigned i = 0;
i < P_order;
i++)
7587 for (
unsigned j = 0; j < P_order; j++)
7589 for (
unsigned k = 0; k < P_order; k++)
7592 dpsids(index, 0) = psi3[
i] * psi2[j] * dpsi1ds[k];
7593 dpsids(index, 1) = psi3[
i] * dpsi2ds[j] * psi1[k];
7594 dpsids(index, 2) = dpsi3ds[
i] * psi2[j] * psi1[k];
7595 psi[index] = psi3[
i] * psi2[j] * psi1[k];
7604 std::ostringstream error_message;
7605 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
7606 error_message <<
"\n polynomial order for shape functions.\n";
7608 OOMPH_CURRENT_FUNCTION,
7609 OOMPH_EXCEPTION_LOCATION);
7617 template<
unsigned INITIAL_NNODE_1D>
7621 std::ostringstream error_message;
7623 <<
"\nd2shape_local currently not implemented for this element\n";
7625 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
7632 template<
unsigned INITIAL_NNODE_1D>
7637 unsigned n_sons = this->tree_pt()->nsons();
7639 unsigned max_son_p_order = 0;
7640 for (
unsigned ison = 0; ison < n_sons; ison++)
7643 this->tree_pt()->son_pt(ison)->object_pt());
7644 son_p_order[ison] = el_pt->
p_order();
7645 if (son_p_order[ison] > max_son_p_order)
7646 max_son_p_order = son_p_order[ison];
7649 unsigned old_Nnode = this->nnode();
7650 unsigned old_P_order = this->p_order();
7652 this->p_order() = max_son_p_order;
7655 delete this->integral_pt();
7656 switch (this->p_order())
7677 std::ostringstream error_message;
7678 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
7679 error_message <<
"\n integration scheme.\n";
7681 error_message.str(),
7682 "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
7683 OOMPH_EXCEPTION_LOCATION);
7688 for (
unsigned n = 0; n < old_Nnode; n++)
7690 old_node_pt[n] = this->node_pt(n);
7694 this->set_n_node(this->p_order() * this->p_order() * this->p_order());
7697 this->node_pt(0) = old_node_pt[0];
7698 this->node_pt(this->p_order() - 1) = old_node_pt[old_P_order - 1];
7699 this->node_pt(this->p_order() * (this->p_order() - 1)) =
7700 old_node_pt[(old_P_order) * (old_P_order - 1)];
7701 this->node_pt(this->p_order() * this->p_order() - 1) =
7702 old_node_pt[(old_P_order) * (old_P_order)-1];
7703 this->node_pt(this->p_order() * this->p_order() * (this->p_order() - 1)) =
7704 old_node_pt[old_P_order * old_P_order * (old_P_order - 1)];
7705 this->node_pt((this->p_order() * this->p_order() + 1) *
7706 (this->p_order() - 1)) =
7707 old_node_pt[(old_P_order * old_P_order + 1) * (old_P_order - 1)];
7708 this->node_pt(this->p_order() * (this->p_order() + 1) *
7709 (this->p_order() - 1)) =
7710 old_node_pt[old_P_order * (old_P_order + 1) * (old_P_order - 1)];
7711 this->node_pt(this->p_order() * this->p_order() * this->p_order() - 1) =
7712 old_node_pt[old_P_order * old_P_order * old_P_order - 1];
7715 if (this->p_order() % 2 == 1)
7718 unsigned n_p = this->p_order();
7719 unsigned m = (n_p - 1) / 2;
7721 unsigned s0space = m;
7722 unsigned s1space = m * n_p;
7723 unsigned s2space = m * n_p * n_p;
7727 this->node_pt(1 * s0space + 0 * s1space + 0 * s2space) =
7730 ->vertex_node_pt(0);
7732 this->node_pt(0 * s0space + 1 * s1space + 0 * s2space) =
7735 ->vertex_node_pt(0);
7737 this->node_pt(1 * s0space + 1 * s1space + 0 * s2space) =
7740 ->vertex_node_pt(0);
7742 this->node_pt(2 * s0space + 1 * s1space + 0 * s2space) =
7745 ->vertex_node_pt(1);
7747 this->node_pt(1 * s0space + 2 * s1space + 0 * s2space) =
7750 ->vertex_node_pt(2);
7754 this->node_pt(0 * s0space + 0 * s1space + 1 * s2space) =
7757 ->vertex_node_pt(0);
7759 this->node_pt(1 * s0space + 0 * s1space + 1 * s2space) =
7762 ->vertex_node_pt(1);
7764 this->node_pt(2 * s0space + 0 * s1space + 1 * s2space) =
7767 ->vertex_node_pt(1);
7769 this->node_pt(0 * s0space + 1 * s1space + 1 * s2space) =
7772 ->vertex_node_pt(0);
7774 this->node_pt(1 * s0space + 1 * s1space + 1 * s2space) =
7777 ->vertex_node_pt(0);
7779 this->node_pt(2 * s0space + 1 * s1space + 1 * s2space) =
7782 ->vertex_node_pt(1);
7784 this->node_pt(0 * s0space + 2 * s1space + 1 * s2space) =
7787 ->vertex_node_pt(2);
7789 this->node_pt(1 * s0space + 2 * s1space + 1 * s2space) =
7792 ->vertex_node_pt(2);
7794 this->node_pt(2 * s0space + 2 * s1space + 1 * s2space) =
7797 ->vertex_node_pt(3);
7801 this->node_pt(1 * s0space + 0 * s1space + 2 * s2space) =
7804 ->vertex_node_pt(5);
7806 this->node_pt(0 * s0space + 1 * s1space + 2 * s2space) =
7809 ->vertex_node_pt(4);
7811 this->node_pt(1 * s0space + 1 * s1space + 2 * s2space) =
7814 ->vertex_node_pt(4);
7816 this->node_pt(2 * s0space + 1 * s1space + 2 * s2space) =
7819 ->vertex_node_pt(5);
7821 this->node_pt(1 * s0space + 2 * s1space + 2 * s2space) =
7824 ->vertex_node_pt(6);
7829 if (this->node_pt(0) == 0)
7832 "The Corner node (0) does not exist",
7833 "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
7834 OOMPH_EXCEPTION_LOCATION);
7837 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
7838 unsigned ntstorage = time_stepper_pt->
ntstorage();
7843 unsigned n_p = this->nnode_1d();
7844 for (
unsigned i0 = 0; i0 < n_p; i0++)
7847 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
7849 s[0] = -1.0 + 2.0 * s_fraction[0];
7851 for (
unsigned i1 = 0; i1 < n_p; i1++)
7854 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
7856 s[1] = -1.0 + 2.0 * s_fraction[1];
7858 for (
unsigned i2 = 0; i2 < n_p; i2++)
7861 s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
7863 s[2] = -1.0 + 2.0 * s_fraction[2];
7866 jnod = i0 + n_p * i1 + n_p * n_p * i2;
7870 bool node_done =
false;
7874 Node* created_node_pt = this->get_node_at_local_coordinate(
s);
7878 if (created_node_pt != 0)
7881 this->node_pt(jnod) = created_node_pt;
7895 bool is_periodic =
false;
7897 node_created_by_neighbour(s_fraction, is_periodic);
7900 if (created_node_pt != 0)
7906 OOMPH_CURRENT_FUNCTION,
7907 OOMPH_EXCEPTION_LOCATION);
7912 this->node_pt(jnod) = created_node_pt;
7929 if (s_fraction[0] < 0.5)
7932 if (s_fraction[1] < 0.5)
7935 if (s_fraction[2] < 0.5)
7939 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
7940 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
7941 s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
7948 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
7949 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
7950 s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
7961 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
7962 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
7963 s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
7970 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
7971 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
7972 s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
7987 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
7988 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
7989 s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
7996 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
7997 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
7998 s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
8009 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
8010 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
8011 s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
8018 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
8019 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
8020 s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
8029 this->tree_pt()->son_pt(son)->object_pt());
8041 else if (i0 == n_p - 1)
8065 else if (i1 == n_p - 1)
8119 else if (i2 == n_p - 1)
8155 std::set<unsigned> boundaries;
8169 if (boundaries.size() > 0)
8173 this->construct_boundary_node(jnod, time_stepper_pt);
8176 Vector<int> bound_cons(this->ncont_interpolated_values());
8177 son_el_pt->
get_bcs(boundary, bound_cons);
8180 unsigned nval = created_node_pt->
nvalue();
8181 for (
unsigned k = 0; k < nval; k++)
8185 created_node_pt->
pin(k);
8192 dynamic_cast<SolidNode*
>(created_node_pt);
8193 if (solid_node_pt != 0)
8196 unsigned n_dim = created_node_pt->
ndim();
8201 if (son_solid_el_pt == 0)
8204 "We have a SolidNode outside a refineable SolidElement\n";
8206 "during mesh refinement -- this doesn't make sense\n";
8209 "PRefineableQElement<3,INITIAL_NNODE_1D>:"
8210 ":rebuild_from_sons()",
8211 OOMPH_EXCEPTION_LOCATION);
8217 for (
unsigned k = 0; k < n_dim; k++)
8219 if (solid_bound_cons[k])
8229 for (std::set<unsigned>::iterator it = boundaries.begin();
8230 it != boundaries.end();
8244 *it, boundary, s_in_son, zeta);
8255 created_node_pt = this->construct_node(jnod, time_stepper_pt);
8272 for (
unsigned t = 0;
t < ntstorage;
t++)
8274 using namespace QuadTreeNames;
8279 son_el_pt->
get_x(
t, s_in_son, x_prev);
8280 for (
unsigned i = 0;
i < 3;
i++)
8282 created_node_pt->
x(
t,
i) = x_prev[
i];
8288 for (
unsigned t = 0;
t < ntstorage;
t++)
8296 for (
unsigned k = 0; k < created_node_pt->
nvalue(); k++)
8298 created_node_pt->
set_value(
t, k, prev_values[k]);
8313 "Have not implemented rebuilding from sons for";
8314 error_message +=
"Algebraic p-refineable elements yet\n";
8318 "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
8319 OOMPH_EXCEPTION_LOCATION);
8335 template<
unsigned INITIAL_NNODE_1D>
8350 template<
unsigned INITIAL_NNODE_1D>
8352 const int& value_id,
const int& my_face, std::ofstream& output_hangfile)
8354 using namespace OcTreeNames;
8359 int neigh_face, diff_level;
8360 bool in_neighbouring_tree;
8370 in_neighbouring_tree);
8377 bool h_type_dependent =
false;
8379 bool p_type_dependent =
false;
8385 if (diff_level != 0)
8388 h_type_dependent =
true;
8391 else if (neigh_pt->
nsons() == 0)
8395 unsigned my_p_order =
8398 unsigned neigh_p_order =
8404 if (neigh_p_order == my_p_order)
8408 else if (neigh_p_order < my_p_order)
8411 p_type_dependent =
true;
8433 if (h_type_dependent || p_type_dependent ||
true)
8441 face_edge_other_face[0] =
B;
8443 face_edge_other_face[1] =
F;
8445 face_edge_other_face[2] =
D;
8447 face_edge_other_face[3] =
U;
8451 face_edge_other_face[0] =
B;
8453 face_edge_other_face[1] =
F;
8455 face_edge_other_face[2] =
D;
8457 face_edge_other_face[3] =
U;
8461 face_edge_other_face[0] =
B;
8463 face_edge_other_face[1] =
F;
8465 face_edge_other_face[2] =
L;
8467 face_edge_other_face[3] =
R;
8471 face_edge_other_face[0] =
B;
8473 face_edge_other_face[1] =
F;
8475 face_edge_other_face[2] =
L;
8477 face_edge_other_face[3] =
R;
8481 face_edge_other_face[0] =
D;
8483 face_edge_other_face[1] =
U;
8485 face_edge_other_face[2] =
L;
8487 face_edge_other_face[3] =
R;
8491 face_edge_other_face[0] =
D;
8493 face_edge_other_face[1] =
U;
8495 face_edge_other_face[2] =
L;
8497 face_edge_other_face[3] =
R;
8501 "my_face not L, R, D, U, B, F\n",
8502 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8503 OOMPH_EXCEPTION_LOCATION);
8507 for (
unsigned i = 0;
i < 4;
i++)
8510 unsigned my_edge = face_edge[
i];
8513 OcTree* edge_neigh_pt = 0;
8517 int neigh_edge = 0, edge_diff_level = 0;
8518 unsigned edge_p_order =
8526 int tmp_neigh_edge, tmp_edge_diff_level;
8529 unsigned i_root_edge_neighbour = 0;
8532 unsigned nroot_edge_neighbour = 0;
8535 bool my_edge_is_master =
true;
8540 bool keep_searching =
true;
8541 bool search_faces =
false;
8542 bool first_face_searched =
false;
8544 while (keep_searching)
8547 OcTree* tmp_edge_neigh_pt;
8553 if (!first_face_searched)
8558 tmp_edge_translate_s,
8559 tmp_edge_s_lo_neigh,
8560 tmp_edge_s_hi_neigh,
8562 tmp_edge_diff_level,
8563 in_neighbouring_tree);
8566 first_face_searched =
true;
8573 tmp_edge_translate_s,
8574 tmp_edge_s_lo_neigh,
8575 tmp_edge_s_hi_neigh,
8577 tmp_edge_diff_level,
8578 in_neighbouring_tree);
8581 keep_searching =
false;
8590 i_root_edge_neighbour,
8591 nroot_edge_neighbour,
8592 tmp_edge_translate_s,
8593 tmp_edge_s_lo_neigh,
8594 tmp_edge_s_hi_neigh,
8596 tmp_edge_diff_level);
8606 bool new_edge_master =
false;
8609 if (tmp_edge_neigh_pt != 0)
8612 if (tmp_edge_diff_level < edge_diff_level)
8616 new_edge_master =
true;
8618 edge_diff_level = tmp_edge_diff_level;
8625 else if (tmp_edge_diff_level == edge_diff_level &&
8626 tmp_edge_neigh_pt->
nsons() == 0)
8633 unsigned tmp_edge_neigh_p_order =
8639 if (tmp_edge_neigh_p_order == edge_p_order)
8643 else if (tmp_edge_neigh_p_order < edge_p_order)
8647 new_edge_master =
true;
8649 edge_diff_level = tmp_edge_diff_level;
8674 if (new_edge_master)
8677 edge_neigh_pt = tmp_edge_neigh_pt;
8678 edge_translate_s = tmp_edge_translate_s;
8679 edge_s_lo_neigh = tmp_edge_s_lo_neigh;
8680 edge_s_hi_neigh = tmp_edge_s_hi_neigh;
8681 neigh_edge = tmp_neigh_edge;
8682 edge_diff_level = tmp_edge_diff_level;
8685 my_edge_is_master =
false;
8690 i_root_edge_neighbour++;
8697 if (i_root_edge_neighbour >= nroot_edge_neighbour)
8701 search_faces =
true;
8707 if (!my_edge_is_master)
8710 unsigned active_coord_index;
8717 active_coord_index = 0;
8723 active_coord_index = 1;
8729 active_coord_index = 2;
8733 "Cannot transform coordinates",
8734 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8735 OOMPH_EXCEPTION_LOCATION);
8751 const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
8752 const unsigned neigh_n_p =
8757 unsigned neighbour_node_number = 0;
8758 Node* neighbour_node_pt = 0;
8761 bool master_is_not_edge =
false;
8762 for (
unsigned i0 = 0; i0 < neigh_n_p; i0++)
8764 const unsigned s0space = 1;
8765 const unsigned s1space = neigh_n_p;
8766 const unsigned s2space = neigh_n_p * neigh_n_p;
8772 neighbour_node_number = i0 * s0space;
8774 neighbour_node_number, value_id);
8777 neighbour_node_number =
8778 (neigh_n_p - 1) * s1space + i0 * s0space;
8780 neighbour_node_number, value_id);
8783 neighbour_node_number =
8784 (neigh_n_p - 1) * s2space + i0 * s0space;
8786 neighbour_node_number, value_id);
8789 neighbour_node_number = (neigh_n_p - 1) * s1space +
8790 (neigh_n_p - 1) * s2space +
8793 neighbour_node_number, value_id);
8797 neighbour_node_number = i0 * s1space;
8799 neighbour_node_number, value_id);
8802 neighbour_node_number =
8803 (neigh_n_p - 1) * s0space + i0 * s1space;
8805 neighbour_node_number, value_id);
8808 neighbour_node_number =
8809 (neigh_n_p - 1) * s2space + i0 * s1space;
8811 neighbour_node_number, value_id);
8814 neighbour_node_number = (neigh_n_p - 1) * s0space +
8815 (neigh_n_p - 1) * s2space +
8818 neighbour_node_number, value_id);
8822 neighbour_node_number = i0 * s2space;
8824 neighbour_node_number, value_id);
8827 neighbour_node_number =
8828 (neigh_n_p - 1) * s0space + i0 * s2space;
8830 neighbour_node_number, value_id);
8833 neighbour_node_number =
8834 (neigh_n_p - 1) * s1space + i0 * s2space;
8836 neighbour_node_number, value_id);
8839 neighbour_node_number = (neigh_n_p - 1) * s0space +
8840 (neigh_n_p - 1) * s1space +
8843 neighbour_node_number, value_id);
8849 master_is_not_edge =
true;
8852 if (master_is_not_edge)
break;
8855 master_node_number.push_back(neighbour_node_number);
8856 master_node_pt.push_back(neighbour_node_pt);
8862 if (master_is_not_edge)
8865 for (
unsigned i0 = 0; i0 < neigh_n_p; i0++)
8868 for (
unsigned i1 = 0; i1 < neigh_n_p; i1++)
8870 const unsigned s0space = 1;
8871 const unsigned s1space = neigh_n_p;
8872 const unsigned s2space = neigh_n_p * neigh_n_p;
8878 neighbour_node_number = i0 * s0space + i1 * s1space;
8881 neighbour_node_number, value_id);
8884 neighbour_node_number =
8885 (neigh_n_p - 1) * s2space + i0 * s0space + i1 * s1space;
8888 neighbour_node_number, value_id);
8891 neighbour_node_number = i0 * s0space + i1 * s2space;
8894 neighbour_node_number, value_id);
8897 neighbour_node_number =
8898 (neigh_n_p - 1) * s1space + i0 * s0space + i1 * s2space;
8901 neighbour_node_number, value_id);
8904 neighbour_node_number = i0 * s1space + i1 * s2space;
8907 neighbour_node_number, value_id);
8910 neighbour_node_number =
8911 (neigh_n_p - 1) * s0space + i0 * s1space + i1 * s2space;
8914 neighbour_node_number, value_id);
8919 "PRefineableQElement<3,INITIAL_NNODE_"
8920 "1D>::oc_hang_helper()",
8921 OOMPH_EXCEPTION_LOCATION);
8925 master_node_number.push_back(neighbour_node_number);
8926 master_node_pt.push_back(neighbour_node_pt);
8933 unsigned local_node_number = 0;
8934 Node* local_node_pt = 0;
8937 for (
unsigned i0 = 0; i0 < my_n_p; i0++)
8942 const unsigned s0space = 1;
8943 const unsigned s1space = my_n_p;
8944 const unsigned s2space = my_n_p * my_n_p;
8952 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
8953 s_fraction[1] = 0.0;
8954 s_fraction[2] = 0.0;
8955 local_node_number = i0 * s0space;
8957 this->interpolating_node_pt(local_node_number, value_id);
8961 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
8962 s_fraction[1] = 1.0;
8963 s_fraction[2] = 0.0;
8964 local_node_number = (my_n_p - 1) * s1space + i0 * s0space;
8966 this->interpolating_node_pt(local_node_number, value_id);
8970 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
8971 s_fraction[1] = 0.0;
8972 s_fraction[2] = 1.0;
8973 local_node_number = (my_n_p - 1) * s2space + i0 * s0space;
8975 this->interpolating_node_pt(local_node_number, value_id);
8979 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
8980 s_fraction[1] = 1.0;
8981 s_fraction[2] = 1.0;
8982 local_node_number = (my_n_p - 1) * s1space +
8983 (my_n_p - 1) * s2space + i0 * s0space;
8985 this->interpolating_node_pt(local_node_number, value_id);
8989 s_fraction[0] = 0.0;
8991 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
8992 s_fraction[2] = 0.0;
8993 local_node_number = i0 * s1space;
8995 this->interpolating_node_pt(local_node_number, value_id);
8998 s_fraction[0] = 1.0;
9000 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9001 s_fraction[2] = 0.0;
9002 local_node_number = (my_n_p - 1) * s0space + i0 * s1space;
9004 this->interpolating_node_pt(local_node_number, value_id);
9007 s_fraction[0] = 0.0;
9009 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9010 s_fraction[2] = 1.0;
9011 local_node_number = (my_n_p - 1) * s2space + i0 * s1space;
9013 this->interpolating_node_pt(local_node_number, value_id);
9016 s_fraction[0] = 1.0;
9018 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9019 s_fraction[2] = 1.0;
9020 local_node_number = (my_n_p - 1) * s0space +
9021 (my_n_p - 1) * s2space + i0 * s1space;
9023 this->interpolating_node_pt(local_node_number, value_id);
9027 s_fraction[0] = 0.0;
9028 s_fraction[1] = 0.0;
9030 local_one_d_fraction_of_interpolating_node(i0, 2, value_id);
9031 local_node_number = i0 * s2space;
9033 this->interpolating_node_pt(local_node_number, value_id);
9036 s_fraction[0] = 1.0;
9037 s_fraction[1] = 0.0;
9039 local_one_d_fraction_of_interpolating_node(i0, 2, value_id);
9040 local_node_number = (my_n_p - 1) * s0space + i0 * s2space;
9042 this->interpolating_node_pt(local_node_number, value_id);
9045 s_fraction[0] = 0.0;
9046 s_fraction[1] = 1.0;
9048 local_one_d_fraction_of_interpolating_node(i0, 2, value_id);
9049 local_node_number = (my_n_p - 1) * s1space + i0 * s2space;
9051 this->interpolating_node_pt(local_node_number, value_id);
9054 s_fraction[0] = 1.0;
9055 s_fraction[1] = 1.0;
9057 local_one_d_fraction_of_interpolating_node(i0, 2, value_id);
9058 local_node_number = (my_n_p - 1) * s0space +
9059 (my_n_p - 1) * s1space + i0 * s2space;
9061 this->interpolating_node_pt(local_node_number, value_id);
9066 "my_edge not recognised\n",
9067 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9068 OOMPH_EXCEPTION_LOCATION);
9072 dependent_node_number.push_back(local_node_number);
9073 dependent_node_pt.push_back(local_node_pt);
9076 dependent_node_s_fraction.push_back(s_fraction);
9080 const unsigned n_dependent_nodes = dependent_node_pt.size();
9081 const unsigned n_master_nodes = master_node_pt.size();
9082 const unsigned dependent_element_nnode_1d = my_n_p;
9083 const unsigned master_element_nnode_1d = neigh_n_p;
9093 dependent_nodal_position,
9098 master_element_nnode_1d, master_nodal_position, master_weight);
9111 const unsigned n_mortar_vertices = 2;
9114 vertex_pos[1] = this->ninterpolating_node_1d(value_id) - 1;
9116 for (
unsigned i = 0;
i < my_n_p - n_mortar_vertices;
i++)
9118 non_vertex_pos[
i] =
i + 1;
9125 std::vector<bool> mortar_vertex_on_master_edge(n_mortar_vertices);
9126 for (
unsigned v = 0; v < n_mortar_vertices; v++)
9129 mortar_vertex_on_master_edge[v] =
true;
9130 unsigned non_extreme_coordinate = 0;
9132 for (
unsigned i = 0;
i < 3;
i++)
9136 edge_s_lo_neigh[
i] +
9137 dependent_node_s_fraction[vertex_pos[v]][edge_translate_s[
i]] *
9138 (edge_s_hi_neigh[
i] - edge_s_lo_neigh[
i]);
9142 if (std::fabs(std::fabs(s_in_neigh[
i]) - 1.0) > 1.0e-14)
9144 non_extreme_coordinate++;
9145 if (non_extreme_coordinate > 1)
9147 mortar_vertex_on_master_edge[v] =
false;
9155 bool my_edge_coincides_with_master =
true;
9156 for (
unsigned v = 0; v < n_mortar_vertices; v++)
9158 my_edge_coincides_with_master =
9159 my_edge_coincides_with_master && mortar_vertex_on_master_edge[v];
9165 for (
unsigned v = 0; v < n_mortar_vertices; v++)
9170 if ((!my_edge_coincides_with_master) &&
9171 mortar_vertex_on_master_edge[v])
9175 bool node_is_shared =
false;
9176 for (
unsigned i = 0;
i < master_node_pt.size();
i++)
9178 if (dependent_node_pt[vertex_pos[v]] == master_node_pt[
i])
9180 node_is_shared =
true;
9187 if (!node_is_shared)
9194 for (
unsigned i = 0;
i < 3;
i++)
9196 s_in_neigh[
i] = edge_s_lo_neigh[
i] +
9197 dependent_node_s_fraction[vertex_pos[v]]
9198 [edge_translate_s[
i]] *
9199 (edge_s_hi_neigh[
i] - edge_s_lo_neigh[
i]);
9204 s_in_neigh, master_shapes, value_id);
9208 dependent_node_pt[vertex_pos[v]]->set_nonhanging();
9212 for (
unsigned m = 0; m < n_master_nodes; m++)
9215 if (std::fabs(master_shapes[master_node_number[m]]) < 1.0e-14)
9218 master_node_zero_weight.push_back(m);
9224 new HangInfo(n_master_nodes - master_node_zero_weight.size());
9225 unsigned mindex = 0;
9226 for (
unsigned m = 0; m < n_master_nodes; m++)
9230 for (
unsigned i = 0;
i < master_node_zero_weight.size();
i++)
9232 if (m == master_node_zero_weight[
i]) skip =
true;
9240 master_shapes[master_node_number[m]]);
9251 dependent_node_pt[vertex_pos[v]]->set_hanging_pt(explicit_hang_pt,
9272 for (
unsigned m = 0; m < explicit_hang_pt->
nmaster(); m++)
9274 if (std::fabs(explicit_hang_pt->
master_weight(m)) < 1.0e-14)
9277 "Master has zero weight!",
9278 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9279 OOMPH_EXCEPTION_LOCATION);
9288 if (n_dependent_nodes - n_mortar_vertices > 0)
9295 for (
unsigned i = 0;
i < shared_node_M.size();
i++)
9297 shared_node_M[
i].resize(n_dependent_nodes - n_mortar_vertices);
9301 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
9305 pow(-1.0,
int((dependent_element_nnode_1d - 1) -
i - 1)) *
9307 dependent_element_nnode_1d - 1,
9308 dependent_nodal_position[non_vertex_pos[
i]]);
9310 diag_M[
i] = psi[
i] * dependent_weight[non_vertex_pos[
i]];
9314 for (
unsigned v = 0; v < shared_node_M.size(); v++)
9316 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
9320 if (std::fabs(dependent_nodal_position[non_vertex_pos[
i]] -
9321 dependent_nodal_position[vertex_pos[v]]) >=
9326 pow(-1.0,
int((dependent_element_nnode_1d - 1) -
i - 1)) *
9328 dependent_element_nnode_1d - 1,
9329 dependent_nodal_position[vertex_pos[v]]) /
9330 (dependent_nodal_position[non_vertex_pos[
i]] -
9331 dependent_nodal_position[vertex_pos[v]]);
9335 dependent_element_nnode_1d - 1,
9336 dependent_nodal_position[vertex_pos[v]])) < 1.0e-8)
9340 pow(-1.0,
int((dependent_element_nnode_1d - 1) -
i - 1)) *
9342 dependent_element_nnode_1d - 1,
9343 dependent_nodal_position[non_vertex_pos[
i]]);
9349 "Cannot use l'Hopital's rule. Dividing by zero is not "
9351 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9352 OOMPH_EXCEPTION_LOCATION);
9355 shared_node_M[v][
i] = psi[
i] * dependent_weight[vertex_pos[v]];
9365 for (
unsigned i = 0;
i < P.size();
i++)
9367 P[
i].resize(n_master_nodes, 0.0);
9382 if (dependent_element_nnode_1d >= master_element_nnode_1d)
9385 quadrature_knot = &dependent_nodal_position;
9386 quadrature_weight = &dependent_weight;
9391 quadrature_knot = &master_nodal_position;
9392 quadrature_weight = &master_weight;
9396 for (
unsigned q = 0; q < (*quadrature_weight).size(); q++)
9402 s_on_mortar[0] = (*quadrature_knot)[q];
9405 for (
unsigned k = 0; k < n_dependent_nodes - n_mortar_vertices;
9409 if (std::fabs(dependent_nodal_position[non_vertex_pos[k]] -
9410 s_on_mortar[0]) >= 1.0e-08)
9414 pow(-1.0,
int((dependent_element_nnode_1d - 1) - k - 1)) *
9417 (dependent_nodal_position[non_vertex_pos[k]] -
9422 dependent_element_nnode_1d - 1, s_on_mortar[0])) <
9427 pow(-1.0,
int((dependent_element_nnode_1d - 1) - k - 1)) *
9435 "Cannot use l'Hopital's rule. Dividing by zero is not "
9437 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9438 OOMPH_EXCEPTION_LOCATION);
9445 for (
unsigned i = 0;
i < 3;
i++)
9447 s_fraction[
i] = (
i == active_coord_index) ?
9448 0.5 * (s_on_mortar[0] + 1.0) :
9449 dependent_node_s_fraction[vertex_pos[0]][
i];
9454 for (
unsigned i = 0;
i < 3;
i++)
9456 s_in_neigh[
i] = edge_s_lo_neigh[
i] +
9457 s_fraction[edge_translate_s[
i]] *
9458 (edge_s_hi_neigh[
i] - edge_s_lo_neigh[
i]);
9463 s_in_neigh, master_shapes, value_id);
9466 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
9469 for (
unsigned j = 0; j < n_master_nodes; j++)
9471 P[
i][j] += master_shapes[master_node_number[j]] * psi[
i] *
9472 (*quadrature_weight)[q];
9494 for (
unsigned v = 0; v < n_mortar_vertices; v++)
9499 for (
unsigned i = 0;
i < 3;
i++)
9502 (
i == active_coord_index) ?
9503 0.5 * (dependent_nodal_position[vertex_pos[v]] + 1.0) :
9504 dependent_node_s_fraction[vertex_pos[0]][
i];
9509 for (
unsigned i = 0;
i < 3;
i++)
9511 s_in_neigh[
i] = edge_s_lo_neigh[
i] +
9512 s_fraction[edge_translate_s[
i]] *
9513 (edge_s_hi_neigh[
i] - edge_s_lo_neigh[
i]);
9518 s_in_neigh, master_shapes, value_id);
9520 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
9523 for (
unsigned k = 0; k < n_master_nodes; k++)
9526 master_shapes[master_node_number[k]] * shared_node_M[v][
i];
9545 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
9547 for (
unsigned j = 0; j < n_master_nodes; j++)
9549 P[
i][j] /= diag_M[
i];
9568 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
9572 for (
unsigned m = 0; m < n_master_nodes; m++)
9575 if (std::fabs(P[
i][m]) < 1.0e-14)
9578 master_node_zero_weight.push_back(m);
9584 new HangInfo(n_master_nodes - master_node_zero_weight.size());
9585 unsigned mindex = 0;
9586 for (
unsigned m = 0; m < n_master_nodes; m++)
9590 for (
unsigned k = 0; k < master_node_zero_weight.size(); k++)
9592 if (m == master_node_zero_weight[k]) skip =
true;
9597 hang_info_pt[
i]->set_master_node_pt(
9598 mindex++, master_node_pt[m], P[
i][m]);
9622 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
9628 bool node_is_really_shared =
false;
9629 for (
unsigned m = 0; m < hang_info_pt[
i]->nmaster(); m++)
9633 if (hang_info_pt[
i]->master_node_pt(m) ==
9634 dependent_node_pt[non_vertex_pos[
i]])
9636 node_is_really_shared =
true;
9640 if (std::fabs(hang_info_pt[
i]->master_weight(m) - 1.0) >
9644 "node at a mortar vertex",
9645 "PRefineableQElemen<2,INITIAL_NNODE_1D>"
9646 "t::quad_hang_helper()",
9647 OOMPH_EXCEPTION_LOCATION);
9654 if (!node_is_really_shared)
9656 dependent_node_pt[non_vertex_pos[
i]]->set_hanging_pt(
9657 hang_info_pt[
i], -1);
9688 if (h_type_dependent || p_type_dependent)
9692 if (my_face ==
B || my_face ==
F)
9694 active_coord_index[0] = 0;
9695 active_coord_index[1] = 1;
9697 else if (my_face ==
D || my_face ==
U)
9699 active_coord_index[0] = 0;
9700 active_coord_index[1] = 2;
9702 else if (my_face ==
L || my_face ==
R)
9704 active_coord_index[0] = 1;
9705 active_coord_index[1] = 2;
9710 "Cannot transform coordinates",
9711 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9712 OOMPH_EXCEPTION_LOCATION);
9727 const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
9732 unsigned neighbour_node_number = 0;
9733 Node* neighbour_node_pt = 0;
9736 for (
unsigned i0 = 0; i0 < neigh_n_p; i0++)
9738 for (
unsigned i1 = 0; i1 < neigh_n_p; i1++)
9740 const unsigned s0space = 1;
9741 const unsigned s1space = neigh_n_p;
9742 const unsigned s2space = neigh_n_p * neigh_n_p;
9748 neighbour_node_number = i0 * s0space + i1 * s1space;
9750 neighbour_node_number, value_id);
9754 neighbour_node_number =
9755 (neigh_n_p - 1) * s2space + i0 * s0space + i1 * s1space;
9757 neighbour_node_number, value_id);
9761 neighbour_node_number = i0 * s0space + i1 * s2space;
9763 neighbour_node_number, value_id);
9767 neighbour_node_number =
9768 (neigh_n_p - 1) * s1space + i0 * s0space + i1 * s2space;
9770 neighbour_node_number, value_id);
9774 neighbour_node_number = i0 * s1space + i1 * s2space;
9776 neighbour_node_number, value_id);
9780 neighbour_node_number =
9781 (neigh_n_p - 1) * s0space + i0 * s1space + i1 * s2space;
9783 neighbour_node_number, value_id);
9788 "my_face not L, R, D, U, B, F\n",
9789 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9790 OOMPH_EXCEPTION_LOCATION);
9794 master_node_number.push_back(neighbour_node_number);
9795 master_node_pt.push_back(neighbour_node_pt);
9800 unsigned local_node_number = 0;
9801 Node* local_node_pt = 0;
9804 for (
unsigned i0 = 0; i0 < my_n_p; i0++)
9806 for (
unsigned i1 = 0; i1 < my_n_p; i1++)
9811 const unsigned s0space = 1;
9812 const unsigned s1space = my_n_p;
9813 const unsigned s2space = my_n_p * my_n_p;
9821 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
9823 local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
9824 s_fraction[2] = 0.0;
9825 local_node_number = i0 * s0space + i1 * s1space;
9827 this->interpolating_node_pt(local_node_number, value_id);
9832 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
9834 local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
9835 s_fraction[2] = 1.0;
9837 (my_n_p - 1) * s2space + i0 * s0space + i1 * s1space;
9839 this->interpolating_node_pt(local_node_number, value_id);
9844 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
9845 s_fraction[1] = 0.0;
9847 local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
9848 local_node_number = i0 * s0space + i1 * s2space;
9850 this->interpolating_node_pt(local_node_number, value_id);
9855 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
9856 s_fraction[1] = 1.0;
9858 local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
9860 (my_n_p - 1) * s1space + i0 * s0space + i1 * s2space;
9862 this->interpolating_node_pt(local_node_number, value_id);
9866 s_fraction[0] = 0.0;
9868 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9870 local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
9871 local_node_number = i0 * s1space + i1 * s2space;
9873 this->interpolating_node_pt(local_node_number, value_id);
9877 s_fraction[0] = 1.0;
9879 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9881 local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
9883 (my_n_p - 1) * s0space + i0 * s1space + i1 * s2space;
9885 this->interpolating_node_pt(local_node_number, value_id);
9890 "my_face not L, R, D, U, B, F\n",
9891 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9892 OOMPH_EXCEPTION_LOCATION);
9896 dependent_node_number.push_back(local_node_number);
9897 dependent_node_pt.push_back(local_node_pt);
9900 dependent_node_s_fraction.push_back(s_fraction);
9905 const unsigned n_dependent_nodes = dependent_node_pt.size();
9906 const unsigned n_master_nodes = master_node_pt.size();
9907 const unsigned dependent_element_nnode_1d = my_n_p;
9908 const unsigned master_element_nnode_1d = neigh_n_p;
9940 dependent_nodal_position_1d,
9941 dependent_weight_1d);
9945 master_element_nnode_1d, master_nodal_position_1d, master_weight_1d);
9949 dependent_element_nnode_1d * dependent_element_nnode_1d);
9950 for (
unsigned i = 0;
i < dependent_nodal_position.size();
i++)
9952 dependent_nodal_position[
i].resize(2);
9955 dependent_element_nnode_1d);
9957 master_element_nnode_1d);
9958 for (
unsigned i = 0;
i < master_nodal_position.size();
i++)
9960 master_nodal_position[
i].resize(2);
9963 master_element_nnode_1d);
9966 unsigned dependent_index = 0;
9967 for (
unsigned i = 0;
i < dependent_element_nnode_1d;
i++)
9969 for (
unsigned j = 0; j < dependent_element_nnode_1d; j++)
9971 dependent_nodal_position[dependent_index][0] =
9972 dependent_nodal_position_1d[
i];
9973 dependent_nodal_position[dependent_index][1] =
9974 dependent_nodal_position_1d[j];
9975 dependent_weight[dependent_index] =
9976 dependent_weight_1d[
i] * dependent_weight_1d[j];
9980 unsigned master_index = 0;
9981 for (
unsigned i = 0;
i < master_element_nnode_1d;
i++)
9983 for (
unsigned j = 0; j < master_element_nnode_1d; j++)
9985 master_nodal_position[master_index][0] = master_nodal_position_1d[
i];
9986 master_nodal_position[master_index][1] = master_nodal_position_1d[j];
9987 master_weight[master_index] =
9988 master_weight_1d[
i] * master_weight_1d[j];
10011 (dependent_element_nnode_1d - 2));
10012 unsigned nvindex = 0;
10013 for (
unsigned i = 1;
i < dependent_element_nnode_1d - 1;
i++)
10015 for (
unsigned j = 1; j < dependent_element_nnode_1d - 1; j++)
10017 non_vertex_pos[nvindex++] =
i * dependent_element_nnode_1d + j;
10021 for (
unsigned i = 0;
i < n_dependent_nodes;
i++)
10024 bool node_is_vertex =
true;
10025 for (
unsigned j = 0; j < non_vertex_pos.size(); j++)
10027 if (
i == non_vertex_pos[j])
10030 node_is_vertex =
false;
10035 if (node_is_vertex)
10037 vertex_pos.push_back(
i);
10041 const unsigned n_mortar_vertices = vertex_pos.size();
10045 if (n_dependent_nodes - n_mortar_vertices > 0)
10052 for (
unsigned i = 0;
i < shared_node_M.size();
i++)
10054 shared_node_M[
i].resize(n_dependent_nodes - n_mortar_vertices);
10058 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
10064 for (
unsigned dir = 0; dir < 2; dir++)
10067 (dir == 0) ?
i :
i % (dependent_element_nnode_1d - 2);
10070 pow(-1.0,
int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10072 dependent_element_nnode_1d - 1,
10073 dependent_nodal_position[non_vertex_pos[
i]][dir]);
10076 diag_M[
i] = psi[
i] * dependent_weight[non_vertex_pos[
i]];
10088 for (
unsigned v = 0; v < shared_node_M.size(); v++)
10090 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
10096 for (
unsigned dir = 0; dir < 2; dir++)
10099 (dir == 0) ?
i :
i % (dependent_element_nnode_1d - 2);
10101 if (std::fabs(dependent_nodal_position[non_vertex_pos[
i]][dir] -
10102 dependent_nodal_position[vertex_pos[v]][dir]) >=
10108 int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10110 dependent_element_nnode_1d - 1,
10111 dependent_nodal_position[vertex_pos[v]][dir]) /
10112 (dependent_nodal_position[non_vertex_pos[
i]][dir] -
10113 dependent_nodal_position[vertex_pos[v]][dir]);
10117 dependent_element_nnode_1d - 1,
10118 dependent_nodal_position[vertex_pos[v]][dir])) <
10124 int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10126 dependent_element_nnode_1d - 1,
10127 dependent_nodal_position[non_vertex_pos[
i]][dir]);
10133 "Cannot use l'Hopital's rule. Dividing by zero is not "
10135 "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10136 OOMPH_EXCEPTION_LOCATION);
10140 shared_node_M[v][
i] = psi[
i] * dependent_weight[vertex_pos[v]];
10161 for (
unsigned i = 0;
i < P.size();
i++)
10163 P[
i].resize(n_master_nodes, 0.0);
10179 if (dependent_element_nnode_1d >= master_element_nnode_1d)
10182 quadrature_knot = &dependent_nodal_position;
10183 quadrature_weight = &dependent_weight;
10188 quadrature_knot = &master_nodal_position;
10189 quadrature_weight = &master_weight;
10193 for (
unsigned q = 0; q < (*quadrature_weight).size(); q++)
10198 for (
unsigned i = 0;
i < 2;
i++)
10200 s_on_mortar[
i] = (*quadrature_knot)[q][
i];
10204 for (
unsigned k = 0; k < n_dependent_nodes - n_mortar_vertices; k++)
10210 for (
unsigned dir = 0; dir < 2; dir++)
10213 (dir == 0) ? k : k % (dependent_element_nnode_1d - 2);
10215 if (std::fabs(dependent_nodal_position[non_vertex_pos[k]][dir] -
10216 s_on_mortar[dir]) >= 1.0e-08)
10221 int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10223 s_on_mortar[dir]) /
10224 (dependent_nodal_position[non_vertex_pos[k]][dir] -
10229 dependent_element_nnode_1d - 1, s_on_mortar[dir])) <
10235 int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10243 "Cannot use l'Hopital's rule. Dividing by zero is not "
10245 "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10246 OOMPH_EXCEPTION_LOCATION);
10253 for (
unsigned i = 0;
i < 3;
i++)
10255 if (
i == active_coord_index[0])
10257 s_fraction[
i] = 0.5 * (s_on_mortar[0] + 1.0);
10259 else if (
i == active_coord_index[1])
10261 s_fraction[
i] = 0.5 * (s_on_mortar[1] + 1.0);
10265 s_fraction[
i] = dependent_node_s_fraction[vertex_pos[0]][
i];
10271 for (
unsigned i = 0;
i < 3;
i++)
10273 s_in_neigh[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]] *
10274 (s_hi_neigh[
i] - s_lo_neigh[
i]);
10279 s_in_neigh, master_shapes, value_id);
10282 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
10284 for (
unsigned j = 0; j < n_master_nodes; j++)
10286 P[
i][j] += master_shapes[master_node_number[j]] * psi[
i] *
10287 (*quadrature_weight)[q];
10309 for (
unsigned v = 0; v < n_mortar_vertices; v++)
10313 for (
unsigned i = 0;
i < 3;
i++)
10315 if (
i == active_coord_index[0])
10318 0.5 * (dependent_nodal_position[vertex_pos[v]][0] + 1.0);
10320 else if (
i == active_coord_index[1])
10323 0.5 * (dependent_nodal_position[vertex_pos[v]][1] + 1.0);
10327 s_fraction[
i] = dependent_node_s_fraction[vertex_pos[0]][
i];
10333 for (
unsigned i = 0;
i < 3;
i++)
10335 s_in_neigh[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]] *
10336 (s_hi_neigh[
i] - s_lo_neigh[
i]);
10341 s_in_neigh, master_shapes, value_id);
10343 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
10345 for (
unsigned k = 0; k < n_master_nodes; k++)
10348 master_shapes[master_node_number[k]] * shared_node_M[v][
i];
10366 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
10368 for (
unsigned j = 0; j < n_master_nodes; j++)
10370 P[
i][j] /= diag_M[
i];
10388 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
10390 hang_info_pt[
i] =
new HangInfo(n_master_nodes);
10395 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
10397 for (
unsigned j = 0; j < n_master_nodes; j++)
10399 hang_info_pt[
i]->set_master_node_pt(j, master_node_pt[j], P[
i][j]);
10405 for (
unsigned i = 0;
i < n_dependent_nodes - n_mortar_vertices;
i++)
10411 bool node_is_really_shared =
false;
10412 for (
unsigned m = 0; m < hang_info_pt[
i]->nmaster(); m++)
10416 if (hang_info_pt[
i]->master_node_pt(m) ==
10417 dependent_node_pt[non_vertex_pos[
i]])
10419 node_is_really_shared =
true;
10423 if (std::fabs(hang_info_pt[
i]->master_weight(m) - 1.0) > 1.0e-06)
10426 "Something fishy here -- with shared node at a mortar vertex",
10427 "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10428 OOMPH_EXCEPTION_LOCATION);
10435 if (!node_is_really_shared)
10437 dependent_node_pt[non_vertex_pos[
i]]->set_hanging_pt(
10438 hang_info_pt[
i], -1);
10446 for (
unsigned i = 0;
i < n_dependent_nodes;
i++)
10449 if (dependent_node_pt[
i]->is_hanging())
10452 double weight_sum = 0.0;
10453 bool zero_weight =
false;
10454 for (
unsigned m = 0;
10455 m < dependent_node_pt[
i]->hanging_pt()->nmaster();
10458 weight_sum += dependent_node_pt[
i]->hanging_pt()->master_weight(m);
10459 if (std::fabs(dependent_node_pt[
i]->hanging_pt()->master_weight(
10462 zero_weight =
true;
10463 oomph_info <<
"In the hanging scheme for dependent node " <<
i
10464 <<
", master node " << m <<
" has weight "
10465 << dependent_node_pt[
i]->hanging_pt()->master_weight(m)
10466 <<
" < 1.0e-14" << std::endl;
10471 if (std::fabs(weight_sum - 1.0) > 1.0e-08)
10473 oomph_info <<
"Sum of master weights: " << weight_sum << std::endl;
10475 "Weights in hanging scheme do not sum to 1",
10476 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10477 OOMPH_EXCEPTION_LOCATION);
10482 "Zero weights present in hanging schemes",
10483 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10484 OOMPH_EXCEPTION_LOCATION);
10489 for (
unsigned m = 0;
10490 m < dependent_node_pt[
i]->hanging_pt()->nmaster();
10493 if (dependent_node_pt[
i] ==
10494 dependent_node_pt[
i]->hanging_pt()->master_node_pt(m))
10498 "Dependent node has itself as a master!",
10499 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10500 OOMPH_EXCEPTION_LOCATION);
10502 if (dependent_node_pt[
i]
10504 ->master_node_pt(m)
10509 dependent_node_pt[
i]->hanging_pt()->master_node_pt(m);
10510 for (
unsigned mm = 0; mm < new_nod_pt->
hanging_pt()->nmaster();
10513 if (dependent_node_pt[
i] ==
10516 std::cout <<
"++++++++++++++++++++++++++++++++++++++++"
10519 <<
" Dependent node: " << dependent_node_pt[
i]
10520 <<
" = " << dependent_node_pt[
i]->x(0) <<
" "
10521 << dependent_node_pt[
i]->x(1) <<
" "
10522 << dependent_node_pt[
i]->x(2) <<
" " << std::endl;
10524 <<
" Master node: "
10525 << dependent_node_pt[
i]->hanging_pt()->master_node_pt(m)
10527 << dependent_node_pt[
i]->hanging_pt()->master_node_pt(m)->x(
10530 << dependent_node_pt[
i]->hanging_pt()->master_node_pt(m)->x(
10533 << dependent_node_pt[
i]->hanging_pt()->master_node_pt(m)->x(
10535 <<
" " << std::endl;
10536 std::cout <<
"Master's master node: "
10537 << dependent_node_pt[
i]
10539 ->master_node_pt(m)
10541 ->master_node_pt(mm)
10543 << dependent_node_pt[
i]
10545 ->master_node_pt(m)
10547 ->master_node_pt(mm)
10550 << dependent_node_pt[
i]
10552 ->master_node_pt(m)
10554 ->master_node_pt(mm)
10557 << dependent_node_pt[
i]
10559 ->master_node_pt(m)
10561 ->master_node_pt(mm)
10563 <<
" " << std::endl;
10567 std::cout <<
"Hanging node: " << dependent_node_pt[
i]->x(0)
10568 <<
" " << dependent_node_pt[
i]->x(1) <<
" "
10569 << dependent_node_pt[
i]->x(2) <<
" " << std::endl;
10570 for (
unsigned m_tmp = 0;
10571 m_tmp < dependent_node_pt[
i]->hanging_pt()->nmaster();
10575 <<
" m = " << m_tmp <<
": "
10576 << dependent_node_pt[
i]
10578 ->master_node_pt(m_tmp)
10581 << dependent_node_pt[
i]
10583 ->master_node_pt(m_tmp)
10586 << dependent_node_pt[
i]
10588 ->master_node_pt(m_tmp)
10592 << dependent_node_pt[
i]->hanging_pt()->master_weight(
10594 <<
" " << std::endl;
10598 std::cout <<
"Master node " << m
10599 <<
" of Hanging node: " << new_nod_pt->
x(0) <<
" "
10600 << new_nod_pt->
x(1) <<
" " << new_nod_pt->
x(2)
10601 <<
" " << std::endl;
10602 for (
unsigned mm_tmp = 0;
10603 mm_tmp < new_nod_pt->
hanging_pt()->nmaster();
10607 <<
" mm = " << mm_tmp <<
": "
10621 "Dependent node has itself as a master of a master!",
10622 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10623 OOMPH_EXCEPTION_LOCATION);
10633 bool is_master =
false;
10634 for (
unsigned n = 0; n < n_master_nodes; n++)
10636 if (dependent_node_pt[
i] == master_node_pt[n])
10668 for (
unsigned i = 0;
i < n_dependent_nodes;
i++)
10671 if (dependent_node_pt[
i]->is_hanging())
10674 if (value_id == -1)
10678 this->local_coordinate_of_node(dependent_node_number[
i], s_local);
10683 this->interpolated_x(s_local, x_in_neighb);
10687 dependent_node_pt[
i]->x(0) = x_in_neighb[0];
10688 dependent_node_pt[
i]->x(1) = x_in_neighb[1];
10689 dependent_node_pt[
i]->x(2) = x_in_neighb[2];
////////////////////////////////////////////////////////////////////
void setup_algebraic_node_update(Node *&node_pt, const Vector< double > &s_father, FiniteElement *father_el_pt) const
Set up node update info for (newly created) algebraic node: I.e. work out its node update information...
BinaryTree class: Recursively defined, generalised binary tree.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
void pin(const unsigned &i)
Pin the i-th stored variable.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
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).
A general Finite Element class.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
unsigned nnode() const
Return the number of nodes.
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...
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Class for multidimensional Gauss Lobatto Legendre integration rules empty - just establishes template...
A Generalised Element class.
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.
void set_master_node_pt(const unsigned &i, Node *const &master_node_pt, const double &weight)
Set the pointer to the i-th master node and its weight.
////////////////////////////////////////////////////////////////////
Vector< GeomObject * > & geom_object_pt()
Vector of (pointers to) geometric objects involved in node update function.
virtual void set_node_update_info(const Vector< GeomObject * > &geom_object_pt)=0
Set node update information: Pass the vector of (pointers to) the geometric objects that affect the n...
////////////////////////////////////////////////////////////////////
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
bool boundary_coordinate_exists(const unsigned &i) const
Indicate whether the i-th boundary has an intrinsic coordinate.
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.
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
virtual void make_periodic(Node *const &node_pt)
Make the node periodic by copying the values from node_pt. Note that the coordinates will always rema...
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
OcTree class: Recursively defined, generalised octree.
static Vector< int > faces_of_common_edge(const int &edge)
Function that, given an edge, returns the two faces on which it.
OcTree * gteq_face_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_sw, Vector< double > &s_ne, int &face, int &diff_level, bool &in_neighbouring_tree) const
Find (pointer to) ‘greater-or-equal-sized face neighbour’ in given direction (L/R/U/D/F/B)....
OcTree * gteq_true_edge_neighbour(const int &direction, const unsigned &i_root_edge_neighbour, unsigned &nroot_edge_neighbour, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level) const
Find (pointer to) ‘greater-or-equal-sized true edge neighbour’ in the given direction (LB,...
Class that returns the shape functions associated with legendre.
static double nodal_position(const unsigned &n)
static void calculate_nodal_positions()
Static function used to populate the stored positions.
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....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
unsigned & p_order()
Access function to P_order.
p-refineable version of RefineableQElement<1,INITIAL_NNODE_1D>. Generic class definitions
Node * get_node_at_local_coordinate(const Vector< double > &s) const
Return the node at the specified local coordinate.
p-refineable version of RefineableQElement<2,INITIAL_NNODE_1D>.
Node * get_node_at_local_coordinate(const Vector< double > &s) const
Return the node at the specified local coordinate.
p-refineable version of RefineableQElement<3,INITIAL_NNODE_1D>.
Node * get_node_at_local_coordinate(const Vector< double > &s) const
Return the node at the specified local coordinate.
A class that is used to template the p-refineable Q elements by dimension. It's really nothing more t...
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
QuadTree class: Recursively defined, generalised quadtree.
QuadTree * gteq_edge_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level, bool &in_neighbouring_tree) const
Return pointer to greater or equal-sized edge neighbour in specified direction; also provide info reg...
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
virtual unsigned ninterpolating_node_1d(const int &value_id)
Return the number of nodes in a one_d direction that are used to interpolate the value_id-th unknown....
virtual unsigned ninterpolating_node(const int &value_id)
Return the number of nodes that are used to interpolate the value_id-th unknown. Default is to assume...
virtual void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
Return the basis functions that are used to interpolate the value_id-th unknown. By default assume is...
virtual Node * interpolating_node_pt(const unsigned &n, const int &value_id)
In mixed elements, different sets of nodes are used to interpolate different unknowns....
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get all continously interpolated function values in this element as a Vector. Note: Vector sets is ow...
Refineable version of QElement<2,NNODE_1D>.
void get_boundaries(const int &edge, std::set< unsigned > &boundaries) const
Determine set of (mesh) boundaries that the element edge/vertex lives on.
void get_bcs(int bound, Vector< int > &bound_cons) const
Determine Vector of boundary conditions along edge (or on vertex) bound (S/W/N/E/SW/SE/NW/NE): For va...
void interpolated_zeta_on_edge(const unsigned &boundary, const int &edge, const Vector< double > &s, Vector< double > &zeta)
Return the value of the intrinsic boundary coordinate interpolated along the edge (S/W/N/E)
Refineable version of QElement<3,NNODE_1D>.
void interpolated_zeta_on_face(const unsigned &boundary, const int &face, const Vector< double > &s, Vector< double > &zeta)
Return the value of the intrinsic boundary coordinate interpolated along the face.
void get_bcs(int bound, Vector< int > &bound_cons) const
Determine Vector of boundary conditions along the element's boundary (or vertex) bound (S/W/N/E/SW/SE...
void get_boundaries(const int &edge, std::set< unsigned > &boundaries) const
Given an element edge/vertex, return a Vector which contains all the (mesh-)boundary numbers that thi...
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Refineable version of Solid quad elements.
void get_solid_bcs(int bound, Vector< int > &solid_bound_cons) const
Determine vector of solid (positional) boundary conditions along edge (or on vertex) bound (S/W/N/E/S...
Refineable version of Solid brick elements.
void get_solid_bcs(int bound, Vector< int > &solid_bound_cons) const
Determine vector of solid (positional) boundary conditions along edge (or on vertex) bound (S/W/N/E/S...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
void pin_position(const unsigned &i)
Pin the nodal position.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
bool is_neighbour_periodic(const int &direction)
Return whether the neighbour in the particular direction is periodic.
A generalised tree base class that abstracts the common functionality between the quad- and octrees u...
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
unsigned nsons() const
Return number of sons (zero if it's a leaf node)
int son_type() const
Return son type.
static const int OMEGA
Default value for an unassigned neighbour.
Tree * son_pt(const int &son_index) const
Return pointer to the son for a given index. Note that to aid code readability specific enums have be...
TreeRoot *& root_pt()
Return pointer to root of the tree.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
double dlegendre(const unsigned &p, const double &x)
Calculates first derivative of Legendre polynomial of degree p at x using three term recursive formul...
double ddlegendre(const unsigned &p, const double &x)
Calculates second derivative of Legendre polynomial of degree p at x using three term recursive formu...
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...