41 template<
unsigned INITIAL_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>
89 template<
unsigned INITIAL_NNODE_1D>
91 const unsigned&
n1d,
const unsigned&
i)
120 std::ostringstream error_message;
121 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
122 error_message <<
"\n shape functions.\n";
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)
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>
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)
243 P_order = initial_p_order;
276 std::ostringstream error_message;
277 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
278 error_message <<
"\n integration scheme.\n";
290 template<
unsigned INITIAL_NNODE_1D>
326 template<
unsigned INITIAL_NNODE_1D>
338 "Cloned copy must be a PRefineableQElement<1,INITIAL_NNODE_1D>!",
355 err_stream <<
"Clone element has a different p-order from me,"
357 <<
"but it is supposed to be a copy!" << std::endl;
369 err_stream <<
"Clone element has a different number of nodes from me,"
371 <<
"but it is supposed to be a copy!" << std::endl;
378 for (
unsigned n = 0;
n < this->
nnode();
n++)
387 err_stream <<
"The nodes belonging to the clone element are different"
389 <<
"from mine, but it is supposed to be a copy!"
405 throw OomphLibError(
"Can't handle generalised nodal positions (yet).",
443 std::ostringstream error_message;
444 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
445 error_message <<
"\n integration scheme.\n";
469 for (
unsigned n = 0;
n < P_order;
n++)
498 for (
unsigned t = 0;
t < ntstorage;
t++)
514 for (
unsigned k = 0;
k <
n_var;
k++)
550 for (
unsigned t = 0;
t < ntstorage;
t++)
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()",
652 std::string error_message =
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";
670 m_el_pt->set_node_update_info(geom_object_pt);
678 for (
unsigned n = 0;
n < P_order;
n++)
686 for (
unsigned t = 0;
t < ntstorage;
t++)
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++)
841 for (
unsigned i = 0;
i < p_order();
i++)
856 for (
unsigned i = 0;
i < p_order();
i++)
871 for (
unsigned i = 0;
i < p_order();
i++)
886 for (
unsigned i = 0;
i < p_order();
i++)
901 for (
unsigned i = 0;
i < p_order();
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";
935 template<
unsigned INITIAL_NNODE_1D>
945 template<
unsigned INITIAL_NNODE_1D>
950 unsigned n_sons = this->tree_pt()->nsons();
956 this->tree_pt()->son_pt(
ison)->object_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";
1012 this->
node_pt(0) = old_node_pt[0];
1079 using namespace BinaryTreeNames;
1098 this->tree_pt()->son_pt(
son)->object_pt());
1106 std::string error_message =
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";
1132 for (
unsigned t = 0;
t < ntstorage;
t++)
1181 std::string error_message =
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()",
1197 template<
unsigned INITIAL_NNODE_1D>
1209 template<
unsigned INITIAL_NNODE_1D>
1251 std::ostringstream error_message;
1252 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1253 error_message <<
"\n shape functions.\n";
1262 template<
unsigned INITIAL_NNODE_1D>
1271 template<
unsigned INITIAL_NNODE_1D>
1273 const unsigned&
n1d,
const unsigned&
i)
1302 std::ostringstream error_message;
1303 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1304 error_message <<
"\n shape functions.\n";
1316 template<
unsigned INITIAL_NNODE_1D>
1327 for (
unsigned i = 0;
i < 2;
i++)
1335 if (std::fabs(
s[
i] + 1.0) <
tol)
1341 else if (std::fabs(
s[
i] - 1.0) <
tol)
1355 if (std::fabs(z[
n] -
s[
i]) <
tol)
1381 template<
unsigned INITIAL_NNODE_1D>
1385 using namespace QuadTreeNames;
1448 "PRefineableQElement<2,INITIAL_NNODE_1D>::node_"
1449 "created_by_neighbour()",
1467 for (
unsigned i = 0;
i < 2;
i++)
1486 quadtree_pt()->root_pt()->is_neighbour_periodic(
edges[
j]);
1505 template<
unsigned INITIAL_NNODE_1D>
1510 using namespace QuadTreeNames;
1574 for (
unsigned i = 0;
i < 2;
i++)
1642 quadtree_pt()->root_pt()->is_neighbour_periodic(
edges[
j]);
1666 template<
unsigned INITIAL_NNODE_1D>
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)
1709 P_order = initial_p_order;
1742 std::ostringstream error_message;
1743 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1744 error_message <<
"\n integration scheme.\n";
1756 template<
unsigned INITIAL_NNODE_1D>
1760 using namespace QuadTreeNames;
1766 if (Father_bound[
n_p].nrow() == 0)
1768 setup_father_bounds();
1775 int son_type = Tree_pt->
son_type();
1783 std::string error_message =
1784 "Something fishy here: I have no father and yet \n";
1785 error_message +=
"I have no nodes. Who has created me then?!\n";
1848 for (
unsigned i = 0;
i < 2;
i++)
1852 0.5 * (
s_lo[
i] + 1.0) *
1856 0.5 * (
s_hi[
i] + 1.0) *
1866 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1930 for (
unsigned t = 0;
t < ntstorage;
t++)
1946 for (
unsigned k = 0;
k <
n_var;
k++)
1980 m_el_pt->set_node_update_info(geom_object_pt);
1989 template<
unsigned INITIAL_NNODE_1D>
2007 unsigned ngeom_object =
m_el_pt->ngeom_object();
2009 for (
unsigned i = 0;
i < ngeom_object;
i++)
2021 "Cloned copy must be a PRefineableQElement<2,INITIAL_NNODE_1D>!",
2038 err_stream <<
"Clone element has a different p-order from me,"
2040 <<
"but it is supposed to be a copy!" << std::endl;
2052 err_stream <<
"Clone element has a different number of nodes from me,"
2054 <<
"but it is supposed to be a copy!" << std::endl;
2061 for (
unsigned n = 0;
n < this->
nnode();
n++)
2070 err_stream <<
"The nodes belonging to the clone element are different"
2072 <<
"from mine, but it is supposed to be a copy!"
2095 for (
unsigned n = 0;
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()",
2127 throw OomphLibError(
"Can't handle generalised nodal positions (yet).",
2165 std::ostringstream error_message;
2166 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
2167 error_message <<
"\n integration scheme.\n";
2195 for (
unsigned i0 = 0;
i0 < P_order;
i0++)
2202 for (
unsigned i1 = 0;
i1 < P_order;
i1++)
2233 for (
unsigned t = 0;
t < ntstorage;
t++)
2249 for (
unsigned k = 0;
k <
n_var;
k++)
2312 "boundaries.size()!=1 seems a bit strange..\n",
2321 error_stream <<
"Periodic node is not on a boundary...\n"
2337 for (
unsigned t = 0;
t < ntstorage;
t++)
2347 for (
unsigned i = 0;
i < 2;
i++)
2358 for (std::set<unsigned>::iterator
it =
boundaries.begin();
2447 "boundaries.size()!=1 seems a bit strange..\n",
2456 error_stream <<
"Periodic node is not on a boundary...\n"
2472 for (
unsigned t = 0;
t < ntstorage;
t++)
2482 for (
unsigned i = 0;
i < 2;
i++)
2493 for (std::set<unsigned>::iterator
it =
boundaries.begin();
2566 throw OomphLibError(
"boundaries.size()!=1 seems a bit strange..\n",
2612 std::string error_message =
2613 "We have a SolidNode outside a refineable SolidElement\n";
2615 "during mesh refinement -- this doesn't make sense";
2625 for (
unsigned k = 0;
k <
n_dim;
k++)
2640 for (std::set<unsigned>::iterator
it =
boundaries.begin();
2683 for (
unsigned t = 0;
t < ntstorage;
t++)
2694 for (
unsigned i = 0;
i < 2;
i++)
2701 for (
unsigned t = 0;
t < ntstorage;
t++)
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()",
2786 std::string error_message =
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";
2804 m_el_pt->set_node_update_info(geom_object_pt);
2858 for (
unsigned i0 = 0;
i0 < P_order;
i0++)
2865 for (
unsigned i1 = 0;
i1 < P_order;
i1++)
2876 for (
unsigned t = 0;
t < ntstorage;
t++)
2887 for (
unsigned i = 0;
i < 2;
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++)
2937 for (
unsigned i = 0;
i < 3;
i++)
2939 for (
unsigned j = 0;
j < 3;
j++)
2954 for (
unsigned i = 0;
i < 4;
i++)
2956 for (
unsigned j = 0;
j < 4;
j++)
2971 for (
unsigned i = 0;
i < 5;
i++)
2973 for (
unsigned j = 0;
j < 5;
j++)
2988 for (
unsigned i = 0;
i < 6;
i++)
2990 for (
unsigned j = 0;
j < 6;
j++)
3005 for (
unsigned i = 0;
i < 7;
i++)
3007 for (
unsigned j = 0;
j < 7;
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";
3027 template<
unsigned INITIAL_NNODE_1D>
3043 for (
unsigned i = 0;
i < 2;
i++)
3045 for (
unsigned j = 0;
j < 2;
j++)
3067 for (
unsigned i = 0;
i < 3;
i++)
3069 for (
unsigned j = 0;
j < 3;
j++)
3091 for (
unsigned i = 0;
i < 4;
i++)
3093 for (
unsigned j = 0;
j < 4;
j++)
3115 for (
unsigned i = 0;
i < 5;
i++)
3117 for (
unsigned j = 0;
j < 5;
j++)
3139 for (
unsigned i = 0;
i < 6;
i++)
3141 for (
unsigned j = 0;
j < 6;
j++)
3163 for (
unsigned i = 0;
i < 7;
i++)
3165 for (
unsigned j = 0;
j < 7;
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";
3191 template<
unsigned INITIAL_NNODE_1D>
3195 std::ostringstream error_message;
3197 <<
"\nd2shape_local currently not implemented for this element\n";
3206 template<
unsigned INITIAL_NNODE_1D>
3210 using namespace QuadTreeNames;
3213 unsigned n_sons = this->tree_pt()->nsons();
3219 this->tree_pt()->son_pt(
ison)->object_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";
3269 this->
set_n_node(this->p_order() * this->p_order());
3272 this->
node_pt(0) = old_node_pt[0];
3274 this->
node_pt(this->p_order() * (this->p_order() - 1)) =
3276 this->
node_pt(this->p_order() * this->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())
3290 this->
node_pt(midpoint * this->p_order()) =
3292 quadtree_pt()->son_pt(SW)->object_pt())
3297 quadtree_pt()->son_pt(NE)->object_pt())
3300 this->
node_pt((midpoint + 1) * this->p_order() - 1) =
3302 quadtree_pt()->son_pt(NE)->object_pt())
3397 using namespace QuadTreeNames;
3443 this->tree_pt()->son_pt(
son)->object_pt());
3455 else if (
i0 ==
n_p - 1)
3479 else if (
i1 ==
n_p - 1)
3523 for (
unsigned k = 0;
k <
nval;
k++)
3545 std::string error_message =
3546 "We have a SolidNode outside a refineable SolidElement\n";
3548 "during mesh refinement -- this doesn't make sense\n";
3558 for (
unsigned k = 0;
k <
n_dim;
k++)
3570 for (std::set<unsigned>::iterator
it =
boundaries.begin();
3613 for (
unsigned t = 0;
t < ntstorage;
t++)
3615 using namespace QuadTreeNames;
3621 for (
unsigned i = 0;
i < 2;
i++)
3629 for (
unsigned t = 0;
t < ntstorage;
t++)
3653 std::string error_message =
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()",
3677 for (
unsigned j = 0;
j < this->
nnode();
j++)
3688 ->set_node_update_info(
3689 this, s_in_node_update_element, geom_object_pt);
3719 for (
unsigned t = 0;
t < ntstorage;
t++)
3730 for (
unsigned i = 0;
i < 2;
i++)
3747 template<
unsigned INITIAL_NNODE_1D>
3767 using namespace QuadTreeNames;
3814 is_periodic = this->tree_pt()->root_pt()->is_neighbour_periodic(
3902 for (
unsigned i = 0;
i < 2;
i++)
3923 for (
int i = 0;
i < 2;
i++)
3953 neigh_pt->object_pt()->get_interpolated_values(
3958 this->get_interpolated_values(
t,
s, values);
3963 neigh_pt->object_pt()->ncont_interpolated_values();
3974 <<
"erru (S)" <<
err <<
" " <<
ival <<
" "
3995 if (max_error > 1
e-9)
3997 oomph_info <<
"\n#------------------------------------ \n#Max error ";
4000 oomph_info <<
"#------------------------------------ \n " << std::endl;
4050 template<
unsigned INITIAL_NNODE_1D>
4054 using namespace QuadTreeNames;
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()",
4247 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
4256 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
4266 local_one_d_fraction_of_interpolating_node(
i0, 1,
value_id);
4275 local_one_d_fraction_of_interpolating_node(
i0, 1,
value_id);
4341 for (
unsigned i = 0;
i < master_node_pt.
size();
i++)
4359 for (
unsigned i = 0;
i < 2;
i++)
4448 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4449 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4463 for (
unsigned i = 0;
i < P.
size();
i++)
4494 for (
unsigned q = 0; q < (*quadrature_weight).size(); q++)
4531 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4532 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4539 for (
unsigned i = 0;
i < 2;
i++)
4548 for (
unsigned i = 0;
i < 2;
i++)
4564 (*quadrature_weight)[q];
4579 for (
unsigned i = 0;
i < 2;
i++)
4589 for (
unsigned i = 0;
i < 2;
i++)
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()",
4687 for (
unsigned m = 0;
4699 "Weights in hanging scheme do not sum to 1",
4700 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
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()",
4776 template<
unsigned INITIAL_NNODE_1D>
4825 std::ostringstream error_message;
4826 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
4827 error_message <<
"\n shape functions.\n";
4836 template<
unsigned INITIAL_NNODE_1D>
4846 template<
unsigned INITIAL_NNODE_1D>
4848 const unsigned&
n1d,
const unsigned&
i)
4877 std::ostringstream error_message;
4878 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
4879 error_message <<
"\n shape functions.\n";
4890 template<
unsigned INITIAL_NNODE_1D>
4901 for (
unsigned i = 0;
i < 3;
i++)
4909 if (std::fabs(
s[
i] + 1.0) <
tol)
4915 else if (std::fabs(
s[
i] - 1.0) <
tol)
4929 if (std::fabs(z[
n] -
s[
i]) <
tol)
4955 template<
unsigned INITIAL_NNODE_1D>
4959 using namespace OcTreeNames;
4970 edges.push_back(LD);
4974 edges.push_back(LB);
4978 edges.push_back(LU);
4982 edges.push_back(LF);
4991 edges.push_back(RD);
4995 edges.push_back(RB);
4999 edges.push_back(RU);
5003 edges.push_back(RF);
5012 edges.push_back(DB);
5016 edges.push_back(DF);
5025 edges.push_back(UB);
5029 edges.push_back(UF);
5103 for (
unsigned i = 0;
i < 3;
i++)
5122 octree_pt()->root_pt()->is_neighbour_periodic(
faces[
j]);
5202 for (
unsigned i = 0;
i < 3;
i++)
5227 is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
5269 template<
unsigned INITIAL_NNODE_1D>
5274 using namespace OcTreeNames;
5285 edges.push_back(LD);
5289 edges.push_back(LB);
5293 edges.push_back(LU);
5297 edges.push_back(LF);
5306 edges.push_back(RD);
5310 edges.push_back(RB);
5314 edges.push_back(RU);
5318 edges.push_back(RF);
5327 edges.push_back(DB);
5331 edges.push_back(DF);
5340 edges.push_back(UB);
5344 edges.push_back(UF);
5401 for (
unsigned i = 0;
i < 3;
i++)
5526 octree_pt()->root_pt()->is_neighbour_periodic(
faces[
j]);
5590 for (
unsigned i = 0;
i < 3;
i++)
5723 is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
5767 template<
unsigned INITIAL_NNODE_1D>
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)
5812 P_order = initial_p_order;
5817 unsigned new_n_node = P_order * P_order * P_order;
5845 std::ostringstream error_message;
5846 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
5847 error_message <<
"\n integration scheme.\n";
5859 template<
unsigned INITIAL_NNODE_1D>
5863 using namespace OcTreeNames;
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();
5886 std::string error_message =
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()",
6016 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
6017 "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
6069 for (
unsigned t = 0;
t < ntstorage;
t++)
6085 for (
unsigned k = 0;
k <
n_var;
k++)
6103 template<
unsigned INITIAL_NNODE_1D>
6118 "Cloned copy must be a PRefineableQElement<3,INITIAL_NNODE_1D>!",
6135 err_stream <<
"Clone element has a different p-order from me,"
6137 <<
"but it is supposed to be a copy!" << std::endl;
6149 err_stream <<
"Clone element has a different number of nodes from me,"
6151 <<
"but it is supposed to be a copy!" << std::endl;
6158 for (
unsigned n = 0;
n < this->
nnode();
n++)
6167 err_stream <<
"The nodes belonging to the clone element are different"
6169 <<
"from mine, but it is supposed to be a copy!"
6185 throw OomphLibError(
"Can't handle generalised nodal positions (yet).",
6198 this->P_order +=
inc;
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";
6232 this->
set_n_node(this->P_order * this->P_order * this->P_order);
6255 for (
unsigned i0 = 0;
i0 < this->P_order;
i0++)
6262 for (
unsigned i1 = 0;
i1 < this->P_order;
i1++)
6269 for (
unsigned i2 = 0;
i2 < this->P_order;
i2++)
6277 jnod =
i0 + this->P_order *
i1 + this->P_order * this->P_order *
i2;
6301 for (
unsigned t = 0;
t < ntstorage;
t++)
6317 for (
unsigned k = 0;
k <
n_var;
k++)
6360 else if (
i0 == this->P_order - 1)
6385 else if (
i1 == this->P_order - 1)
6440 else if (
i2 == this->P_order - 1)
6494 "boundaries.size()>2 seems a bit strange..\n",
6503 error_stream <<
"Periodic node is not on a boundary...\n"
6519 for (
unsigned t = 0;
t < ntstorage;
t++)
6529 for (
unsigned i = 0;
i <
n_dim;
i++)
6540 for (std::set<unsigned>::iterator
it =
boundaries.begin();
6608 else if (
i0 == this->P_order - 1)
6633 else if (
i1 == this->P_order - 1)
6688 else if (
i2 == this->P_order - 1)
6742 "boundaries.size()>2 seems a bit strange..\n",
6751 error_stream <<
"Periodic node is not on a boundary...\n"
6767 for (
unsigned t = 0;
t < ntstorage;
t++)
6777 for (
unsigned i = 0;
i <
n_dim;
i++)
6788 for (std::set<unsigned>::iterator
it =
boundaries.begin();
6840 else if (
i0 == this->P_order - 1)
6864 else if (
i1 == this->P_order - 1)
6918 else if (
i2 == this->P_order - 1)
6970 "boundaries.size()>2 seems a bit strange..\n",
6971 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
7016 std::string error_message =
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()",
7030 for (
unsigned k = 0;
k <
n_dim;
k++)
7045 for (std::set<unsigned>::iterator
it =
boundaries.begin();
7088 for (
unsigned t = 0;
t < ntstorage;
t++)
7099 for (
unsigned i = 0;
i < 3;
i++)
7106 for (
unsigned t = 0;
t < ntstorage;
t++)
7133 std::string error_message =
7134 "Have not implemented p-refinement for";
7135 error_message +=
"Algebraic p-refineable elements yet\n";
7139 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
7178 for (
unsigned i0 = 0;
i0 < this->P_order;
i0++)
7185 for (
unsigned i1 = 0;
i1 < this->P_order;
i1++)
7192 for (
unsigned i2 = 0;
i2 < this->P_order;
i2++)
7200 jnod =
i0 + this->P_order *
i1 + this->P_order * this->P_order *
i2;
7203 for (
unsigned t = 0;
t < ntstorage;
t++)
7214 for (
unsigned i = 0;
i < 3;
i++)
7245 std::string error_message =
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()",
7265 m_el_pt->set_node_update_info(geom_object_pt);
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) =
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) =
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) =
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) =
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) =
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) =
7411 std::ostringstream error_message;
7412 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
7413 error_message <<
"\n polynomial order for shape functions.\n";
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++)
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++)
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++)
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++)
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++)
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++)
7604 std::ostringstream error_message;
7605 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
7606 error_message <<
"\n polynomial order for shape functions.\n";
7617 template<
unsigned INITIAL_NNODE_1D>
7621 std::ostringstream error_message;
7623 <<
"\nd2shape_local currently not implemented for this element\n";
7632 template<
unsigned INITIAL_NNODE_1D>
7637 unsigned n_sons = this->tree_pt()->nsons();
7643 this->tree_pt()->son_pt(
ison)->object_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()",
7694 this->
set_n_node(this->p_order() * this->p_order() * this->p_order());
7697 this->
node_pt(0) = old_node_pt[0];
7699 this->
node_pt(this->p_order() * (this->p_order() - 1)) =
7701 this->
node_pt(this->p_order() * this->p_order() - 1) =
7703 this->
node_pt(this->p_order() * this->p_order() * (this->p_order() - 1)) =
7705 this->
node_pt((this->p_order() * this->p_order() + 1) *
7706 (this->p_order() - 1)) =
7708 this->
node_pt(this->p_order() * (this->p_order() + 1) *
7709 (this->p_order() - 1)) =
7711 this->
node_pt(this->p_order() * this->p_order() * this->p_order() - 1) =
7715 if (this->p_order() % 2 == 1)
7718 unsigned n_p = this->p_order();
7719 unsigned m = (
n_p - 1) / 2;
7832 "The Corner node (0) does not exist",
7833 "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
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)
8181 for (
unsigned k = 0;
k <
nval;
k++)
8203 std::string error_message =
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()",
8217 for (
unsigned k = 0;
k <
n_dim;
k++)
8229 for (std::set<unsigned>::iterator
it =
boundaries.begin();
8272 for (
unsigned t = 0;
t < ntstorage;
t++)
8274 using namespace QuadTreeNames;
8280 for (
unsigned i = 0;
i < 3;
i++)
8288 for (
unsigned t = 0;
t < ntstorage;
t++)
8312 std::string error_message =
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()",
8335 template<
unsigned INITIAL_NNODE_1D>
8350 template<
unsigned INITIAL_NNODE_1D>
8354 using namespace OcTreeNames;
8501 "my_face not L, R, D, U, B, F\n",
8502 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8507 for (
unsigned i = 0;
i < 4;
i++)
8557 this->octree_pt()->gteq_face_neighbour(
my_face,
8589 this->octree_pt()->gteq_true_edge_neighbour(
my_edge,
8733 "Cannot transform coordinates",
8734 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8919 "PRefineableQElement<3,INITIAL_NNODE_"
8920 "1D>::oc_hang_helper()",
8952 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
8961 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
8970 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
8979 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
8991 local_one_d_fraction_of_interpolating_node(
i0, 1,
value_id);
9000 local_one_d_fraction_of_interpolating_node(
i0, 1,
value_id);
9009 local_one_d_fraction_of_interpolating_node(
i0, 1,
value_id);
9018 local_one_d_fraction_of_interpolating_node(
i0, 1,
value_id);
9030 local_one_d_fraction_of_interpolating_node(
i0, 2,
value_id);
9039 local_one_d_fraction_of_interpolating_node(
i0, 2,
value_id);
9048 local_one_d_fraction_of_interpolating_node(
i0, 2,
value_id);
9057 local_one_d_fraction_of_interpolating_node(
i0, 2,
value_id);
9066 "my_edge not recognised\n",
9067 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9132 for (
unsigned i = 0;
i < 3;
i++)
9142 if (std::fabs(std::fabs(
s_in_neigh[
i]) - 1.0) > 1.0e-14)
9176 for (
unsigned i = 0;
i < master_node_pt.
size();
i++)
9194 for (
unsigned i = 0;
i < 3;
i++)
9277 "Master has zero weight!",
9278 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9349 "Cannot use l'Hopital's rule. Dividing by zero is not "
9351 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9365 for (
unsigned i = 0;
i < P.
size();
i++)
9396 for (
unsigned q = 0; q < (*quadrature_weight).size(); q++)
9435 "Cannot use l'Hopital's rule. Dividing by zero is not "
9437 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9445 for (
unsigned i = 0;
i < 3;
i++)
9454 for (
unsigned i = 0;
i < 3;
i++)
9472 (*quadrature_weight)[q];
9499 for (
unsigned i = 0;
i < 3;
i++)
9509 for (
unsigned i = 0;
i < 3;
i++)
9575 if (std::fabs(P[
i][
m]) < 1.0e-14)
9644 "node at a mortar vertex",
9645 "PRefineableQElemen<2,INITIAL_NNODE_1D>"
9646 "t::quad_hang_helper()",
9710 "Cannot transform coordinates",
9711 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9788 "my_face not L, R, D, U, B, F\n",
9789 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9821 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
9823 local_one_d_fraction_of_interpolating_node(
i1, 1,
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);
9844 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
9847 local_one_d_fraction_of_interpolating_node(
i1, 2,
value_id);
9855 local_one_d_fraction_of_interpolating_node(
i0, 0,
value_id);
9858 local_one_d_fraction_of_interpolating_node(
i1, 2,
value_id);
9868 local_one_d_fraction_of_interpolating_node(
i0, 1,
value_id);
9870 local_one_d_fraction_of_interpolating_node(
i1, 2,
value_id);
9879 local_one_d_fraction_of_interpolating_node(
i0, 1,
value_id);
9881 local_one_d_fraction_of_interpolating_node(
i1, 2,
value_id);
9890 "my_face not L, R, D, U, B, F\n",
9891 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10133 "Cannot use l'Hopital's rule. Dividing by zero is not "
10135 "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10161 for (
unsigned i = 0;
i < P.
size();
i++)
10193 for (
unsigned q = 0; q < (*quadrature_weight).size(); q++)
10198 for (
unsigned i = 0;
i < 2;
i++)
10243 "Cannot use l'Hopital's rule. Dividing by zero is not "
10245 "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10253 for (
unsigned i = 0;
i < 3;
i++)
10271 for (
unsigned i = 0;
i < 3;
i++)
10287 (*quadrature_weight)[q];
10313 for (
unsigned i = 0;
i < 3;
i++)
10333 for (
unsigned i = 0;
i < 3;
i++)
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()",
10454 for (
unsigned m = 0;
10463 oomph_info <<
"In the hanging scheme for dependent node " <<
i
10464 <<
", master node " <<
m <<
" has weight "
10466 <<
" < 1.0e-14" << std::endl;
10475 "Weights in hanging scheme do not sum to 1",
10476 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10482 "Zero weights present in hanging schemes",
10483 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10489 for (
unsigned m = 0;
10498 "Dependent node has itself as a master!",
10499 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10504 ->master_node_pt(
m)
10516 std::cout <<
"++++++++++++++++++++++++++++++++++++++++"
10524 <<
" Master node: "
10535 <<
" " << std::endl;
10536 std::cout <<
"Master's master node: "
10539 ->master_node_pt(
m)
10541 ->master_node_pt(
mm)
10545 ->master_node_pt(
m)
10547 ->master_node_pt(
mm)
10552 ->master_node_pt(
m)
10554 ->master_node_pt(
mm)
10559 ->master_node_pt(
m)
10561 ->master_node_pt(
mm)
10563 <<
" " << std::endl;
10570 for (
unsigned m_tmp = 0;
10575 <<
" m = " <<
m_tmp <<
": "
10578 ->master_node_pt(
m_tmp)
10583 ->master_node_pt(
m_tmp)
10588 ->master_node_pt(
m_tmp)
10594 <<
" " << std::endl;
10598 std::cout <<
"Master node " <<
m
10599 <<
" of Hanging node: " <<
new_nod_pt->x(0) <<
" "
10601 <<
" " << std::endl;
10602 for (
unsigned mm_tmp = 0;
10607 <<
" mm = " <<
mm_tmp <<
": "
10621 "Dependent node has itself as a master of a master!",
10622 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
Base class for algebraic elements.
BinaryTree class: Recursively defined, generalised binary tree.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
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.
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Set pointer to macro element – can be overloaded in derived elements to perform additional tasks.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
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.
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
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.
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
void set_n_node(const unsigned &n)
Set the number of nodes in the element to n, by resizing the storage for pointers to the Node objects...
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
virtual double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
Get the local fraction of any node in the n-th position in a one dimensional expansion along the i-th...
int get_node_number(Node *const &node_pt) const
Return the number of the node *node_pt if this node is in the element, else return -1;.
A Generalised Element class.
unsigned ndim() const
Access function to # of Eulerian coordinates.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Class that contains data for hanging nodes.
Base class for elements that allow MacroElement-based node update.
MacroElementNodeUpdate nodes are nodes with a positional update function, based on their element's Ma...
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.
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.
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....
p-refineable version of RefineableElement
A class that is used to template the p-refineable Q elements by dimension. It's really nothing more t...
QuadTree class: Recursively defined, generalised quadtree.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
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)....
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned nvertex_node() const
Number of vertex nodes in the element.
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
A generalised tree base class that abstracts the common functionality between the quad- and octrees u...
int son_type() const
Return son type.
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
static const int OMEGA
Default value for an unassigned neighbour.
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.
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...