46 template<
unsigned INITIAL_NNODE_1D>
52 switch(this->nnode_1d())
79 oomph_info <<
"\n ERROR: Exceeded maximum polynomial order for";
80 oomph_info <<
"\n shape functions." << std::endl;
86 template<
unsigned INITIAL_NNODE_1D>
90 this->local_coordinate_of_node(n,s_fraction);
91 s_fraction[0] = 0.5*(s_fraction[0] + 1.0);
94 template<
unsigned INITIAL_NNODE_1D>
98 switch(this->nnode_1d())
119 std::ostringstream error_message;
120 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
121 error_message <<
"\n shape functions.\n";
123 OOMPH_CURRENT_FUNCTION,
124 OOMPH_EXCEPTION_LOCATION);
133 template<
unsigned INITIAL_NNODE_1D>
146 if(std::fabs(s[0] + 1.0) < tol)
151 else if(std::fabs(s[0] - 1.0) < tol)
153 index[0] = this->nnode_1d()-1;
162 for (
unsigned n=1; n<this->nnode_1d()-1; n++)
164 if (std::fabs(z[n] - s[0]) < tol)
174 return this->node_pt(index[0]);
184 template<
unsigned INITIAL_NNODE_1D>
199 template<
unsigned INITIAL_NNODE_1D>
207 if (adopted_father_pt!=0)
210 father_pt =
dynamic_cast<BinaryTree*
>(adopted_father_pt);
216 father_pt =
dynamic_cast<BinaryTree*
>(binary_tree_pt()->father_pt());
227 if (father_pt!=0 || initial_p_order!=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>
329 Mesh*
const &mesh_pt,
340 "Cloned copy must be a PRefineableQElement<1,INITIAL_NNODE_1D>!",
341 OOMPH_CURRENT_FUNCTION,
342 OOMPH_EXCEPTION_LOCATION);
349 bool clone_is_ok =
true;
352 clone_is_ok = clone_is_ok && (clone_el_pt->
p_order() == this->p_order());
356 std::ostringstream err_stream;
357 err_stream <<
"Clone element has a different p-order from me," 359 <<
"but it is supposed to be a copy!" << std::endl;
362 OOMPH_CURRENT_FUNCTION,
363 OOMPH_EXCEPTION_LOCATION);
367 clone_is_ok = clone_is_ok && (clone_el_pt->
nnode() == this->nnode());
371 std::ostringstream err_stream;
372 err_stream <<
"Clone element has a different number of nodes from me," 374 <<
"but it is supposed to be a copy!" << std::endl;
377 OOMPH_CURRENT_FUNCTION,
378 OOMPH_EXCEPTION_LOCATION);
382 for(
unsigned n=0; n<this->nnode(); n++)
384 clone_is_ok = clone_is_ok && (clone_el_pt->
node_pt(n) == this->node_pt(n));
389 std::ostringstream err_stream;
390 err_stream <<
"The nodes belonging to the clone element are different" 392 <<
"from mine, but it is supposed to be a copy!" << std::endl;
395 OOMPH_CURRENT_FUNCTION,
396 OOMPH_EXCEPTION_LOCATION);
409 "Can't handle generalised nodal positions (yet).",
410 OOMPH_CURRENT_FUNCTION,
411 OOMPH_EXCEPTION_LOCATION);
416 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
419 unsigned ntstorage = time_stepper_pt->
ntstorage();
425 delete this->integral_pt();
447 std::ostringstream error_message;
448 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
449 error_message <<
"\n integration scheme.\n";
451 OOMPH_CURRENT_FUNCTION,
452 OOMPH_EXCEPTION_LOCATION);
456 this->set_n_node(P_order);
473 for(
unsigned n=0;n<P_order;n++)
477 s_fraction[0] = this->local_one_d_fraction_of_node(n,0);
480 s[0] = s_left[0] + (s_right[0]-s_left[0])*s_fraction[0];
483 bool node_done =
false;
491 if(created_node_pt!=0)
494 this->node_pt(n) = created_node_pt;
502 for(
unsigned t=0;
t<ntstorage;
t++)
511 unsigned n_val_at_node = created_node_pt->
nvalue();
512 unsigned n_val_from_function = prev_values.size();
514 unsigned n_var = n_val_at_node < n_val_from_function ?
515 n_val_at_node : n_val_from_function;
517 for(
unsigned k=0;k<n_var;k++)
519 created_node_pt->
set_value(
t,k,prev_values[k]);
538 created_node_pt = this->construct_node(n,time_stepper_pt);
553 for(
unsigned t=0;
t<ntstorage;
t++)
560 clone_el_pt->
get_x(
t,s,x_prev);
563 created_node_pt->
x(
t,0) = x_prev[0];
577 const unsigned n_value = created_node_pt->
nvalue();
580 for(
unsigned v=0;v<n_value;v++)
582 created_node_pt->
set_value(
t,v,prev_values[v]);
599 "Have not implemented p-refinement for";
601 "Algebraic p-refineable elements yet\n";
605 "PRefineableQElement<1,INITIAL_NNODE_1D>::p_refine()",
606 OOMPH_EXCEPTION_LOCATION);
644 if (clone_m_el_pt!=0)
658 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
660 "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
661 error_message +=
"then I should be too....\n";
664 OOMPH_CURRENT_FUNCTION,
665 OOMPH_EXCEPTION_LOCATION);
676 m_el_pt->set_node_update_info(geom_object_pt);
684 for(
unsigned n=0;n<P_order;n++)
687 s_fraction[0] = this->local_one_d_fraction_of_node(n,0);
689 s[0] = s_left[0] + (s_right[0]-s_left[0])*s_fraction[0];
692 for (
unsigned t=0;
t<ntstorage;
t++)
700 this->get_x(
t,s,x_prev);
703 this->node_pt(n)->x(
t,0) = x_prev[0];
711 this->further_build();
717 template<
unsigned INITIAL_NNODE_1D>
731 for(
unsigned i=0;
i<p_order();
i++) {psi(
i) = psi1[
i];}
742 for(
unsigned i=0;
i<p_order();
i++) {psi(
i) = psi1[
i];}
753 for(
unsigned i=0;
i<p_order();
i++) {psi(
i) = psi1[
i];}
764 for(
unsigned i=0;
i<p_order();
i++) {psi(
i) = psi1[
i];}
775 for(
unsigned i=0;
i<p_order();
i++) {psi(
i) = psi1[
i];}
786 for(
unsigned i=0;
i<p_order();
i++) {psi(
i) = psi1[
i];}
790 oomph_info <<
"\n ERROR: PRefineableQElement::shape() exceeded maximum";
791 oomph_info <<
"\n polynomial order for shape functions." << std::endl;
799 template<
unsigned INITIAL_NNODE_1D>
813 for(
unsigned i=0;
i<p_order();
i++)
816 dpsi(
i,0) = dpsi1ds[
i];
829 for(
unsigned i=0;
i<p_order();
i++)
832 dpsi(
i,0) = dpsi1ds[
i];
844 for(
unsigned i=0;
i<p_order();
i++)
847 dpsi(
i,0) = dpsi1ds[
i];
859 for(
unsigned i=0;
i<p_order();
i++)
862 dpsi(
i,0) = dpsi1ds[
i];
874 for(
unsigned i=0;
i<p_order();
i++)
877 dpsi(
i,0) = dpsi1ds[
i];
889 for(
unsigned i=0;
i<p_order();
i++)
892 dpsi(
i,0) = dpsi1ds[
i];
897 oomph_info <<
"\n ERROR: PRefineableQElement::dshape_local() exceeded maximum";
898 oomph_info <<
"\n polynomial order for shape functions." << std::endl;
907 template<
unsigned INITIAL_NNODE_1D>
912 std::ostringstream error_message;
913 error_message <<
"\nd2shape_local currently not implemented for this element\n";
915 OOMPH_CURRENT_FUNCTION,
916 OOMPH_EXCEPTION_LOCATION);
923 template<
unsigned INITIAL_NNODE_1D>
926 std::ofstream& output_hangfile)
933 template<
unsigned INITIAL_NNODE_1D>
938 unsigned n_sons = this->tree_pt()->nsons();
940 unsigned max_son_p_order = 0;
941 for (
unsigned ison=0;ison<n_sons;ison++)
944 son_p_order[ison] = el_pt->
p_order();
945 if (son_p_order[ison] > max_son_p_order) max_son_p_order = son_p_order[ison];
948 unsigned old_Nnode = this->nnode();
949 unsigned old_P_order = this->p_order();
951 this->p_order() = max_son_p_order;
954 delete this->integral_pt();
955 switch(this->p_order())
976 std::ostringstream error_message;
977 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
978 error_message <<
"\n integration scheme.\n";
980 OOMPH_CURRENT_FUNCTION,
981 OOMPH_EXCEPTION_LOCATION);
986 for (
unsigned n=0; n<old_Nnode; n++)
988 old_node_pt[n] = this->node_pt(n);
992 this->set_n_node(this->p_order());
998 this->node_pt(0) = old_node_pt[0];
999 this->node_pt(this->p_order()-1) = old_node_pt[old_P_order-1];
1007 if(this->node_pt(0)==0)
1010 OOMPH_CURRENT_FUNCTION,
1011 OOMPH_EXCEPTION_LOCATION);
1014 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
1017 const unsigned ntstorage = time_stepper_pt->
ntstorage();
1023 const unsigned n_node = this->nnode_1d();
1026 for(
unsigned n=0;n<n_node;n++)
1029 s_fraction[0] = this->local_one_d_fraction_of_node(n,0);
1032 s[0] = -1.0 + 2.0*s_fraction[0];
1035 if(this->node_pt(n)==0)
1038 bool is_periodic =
false;
1039 Node* created_node_pt =
1040 this->node_created_by_neighbour(s_fraction,is_periodic);
1043 if(created_node_pt!=0)
1049 "Cannot handle periodic nodes yet",
1050 OOMPH_CURRENT_FUNCTION,
1051 OOMPH_EXCEPTION_LOCATION);
1056 this->node_pt(n) = created_node_pt;
1066 using namespace BinaryTreeNames;
1069 if(s_fraction[0] < 0.5)
1072 s_in_son[0] = -1.0 + 4.0*s_fraction[0];
1078 s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
1085 this->tree_pt()->son_pt(son)->object_pt());
1091 if(n==0 || n==n_node-1)
1094 "I am trying to rebuild one of the (two) vertex nodes in\n";
1096 "this 1D element. It should not have been possible to delete\n";
1098 "either of these!\n";
1102 OOMPH_CURRENT_FUNCTION,
1103 OOMPH_EXCEPTION_LOCATION);
1108 this->node_pt(n) = this->construct_node(n,time_stepper_pt);
1122 for(
unsigned t=0;
t<ntstorage;
t++)
1129 son_el_pt->
get_x(
t,s_in_son,x_prev);
1132 this->node_pt(n)->x(
t,0) = x_prev[0];
1146 const unsigned n_value = this->node_pt(n)->nvalue();
1149 for(
unsigned v=0;v<n_value;v++)
1151 this->node_pt(n)->set_value(
t,v,prev_values[v]);
1173 "Have not implemented rebuilding from sons for";
1175 "Algebraic p-refineable elements yet\n";
1179 "PRefineableQElement<1,INITIAL_NNODE_1D>::rebuild_from_sons()",
1180 OOMPH_EXCEPTION_LOCATION);
1190 template<
unsigned INITIAL_NNODE_1D>
1202 template<
unsigned INITIAL_NNODE_1D>
1207 unsigned Nnode_1d = this->nnode_1d();
1208 unsigned n0 = n%Nnode_1d;
1209 unsigned n1 = n/Nnode_1d;
1244 std::ostringstream error_message;
1245 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1246 error_message <<
"\n shape functions.\n";
1248 OOMPH_CURRENT_FUNCTION,
1249 OOMPH_EXCEPTION_LOCATION);
1255 template<
unsigned INITIAL_NNODE_1D>
1259 this->local_coordinate_of_node(n,s_fraction);
1260 s_fraction[0] = 0.5*(s_fraction[0] + 1.0);
1261 s_fraction[1] = 0.5*(s_fraction[1] + 1.0);
1264 template<
unsigned INITIAL_NNODE_1D>
1268 switch(this->nnode_1d())
1289 std::ostringstream error_message;
1290 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1291 error_message <<
"\n shape functions.\n";
1293 OOMPH_CURRENT_FUNCTION,
1294 OOMPH_EXCEPTION_LOCATION);
1303 template<
unsigned INITIAL_NNODE_1D>
1307 unsigned Nnode_1d = this->nnode_1d();
1314 for (
unsigned i=0;
i<2;
i++)
1319 bool is_found=
false;
1322 if(std::fabs(s[
i] + 1.0) < tol)
1328 else if(std::fabs(s[
i] - 1.0) < tol)
1330 index[
i] = Nnode_1d-1;
1340 for (
unsigned n=1; n<Nnode_1d-1; n++)
1342 if (std::fabs(z[n] - s[
i]) < tol)
1358 return this->node_pt(index[0] + Nnode_1d*index[1]);
1368 template<
unsigned INITIAL_NNODE_1D>
1372 using namespace QuadTreeNames;
1376 if(s_fraction[0]==0.0) {edges.push_back(
W);}
1377 if(s_fraction[0]==1.0) {edges.push_back(
E);}
1378 if(s_fraction[1]==0.0) {edges.push_back(
S);}
1379 if(s_fraction[1]==1.0) {edges.push_back(
N);}
1382 unsigned n_size = edges.size();
1384 if(n_size==0) {
return 0;}
1391 int neigh_edge, diff_level;
1392 bool in_neighbouring_tree;
1395 for(
unsigned j=0;j<n_size;j++)
1399 neigh_pt = quadtree_pt()->
1400 gteq_edge_neighbour(edges[j],translate_s,s_lo_neigh,s_hi_neigh,
1401 neigh_edge,diff_level,in_neighbouring_tree);
1410 bool a_vertex_node_is_built =
false;
1417 "PRefineableQElement<2,INITIAL_NNODE_1D>::node_created_by_neighbour()",
1418 OOMPH_EXCEPTION_LOCATION);
1420 for(
unsigned vnode=0; vnode<neigh_obj_pt->nvertex_node(); vnode++)
1422 if(neigh_obj_pt->vertex_node_pt(vnode)!=0)
1423 a_vertex_node_is_built =
true;
1426 if(a_vertex_node_is_built)
1435 for(
unsigned i=0;
i<2;
i++)
1437 s[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
1438 (s_hi_neigh[
i] - s_lo_neigh[
i]);
1442 Node* neighbour_node_pt =
1446 if(neighbour_node_pt!=0)
1450 if(in_neighbouring_tree)
1454 quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1457 return neighbour_node_pt;
1473 template<
unsigned INITIAL_NNODE_1D>
1478 using namespace QuadTreeNames;
1482 if(s_fraction[0]==0.0) {edges.push_back(
W);}
1483 if(s_fraction[0]==1.0) {edges.push_back(
E);}
1484 if(s_fraction[1]==0.0) {edges.push_back(
S);}
1485 if(s_fraction[1]==1.0) {edges.push_back(
N);}
1488 unsigned n_size = edges.size();
1490 if(n_size==0) {
return 0;}
1497 int neigh_edge, diff_level;
1498 bool in_neighbouring_tree;
1501 for(
unsigned j=0;j<n_size;j++)
1505 neigh_pt = quadtree_pt()->
1506 gteq_edge_neighbour(edges[j],translate_s,s_lo_neigh,s_hi_neigh,
1507 neigh_edge,diff_level,in_neighbouring_tree);
1523 for(
unsigned i=0;
i<2;
i++)
1525 s[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
1526 (s_hi_neigh[
i] - s_lo_neigh[
i]);
1530 if(neigh_pt->
nsons()!=0)
1545 s_in_son[0] = 1.0 + 2.0*s[0];
1546 s_in_son[1] = 1.0 + 2.0*s[1];
1553 s_in_son[0] = 1.0 + 2.0*s[0];
1554 s_in_son[1] = -1.0 + 2.0*s[1];
1564 s_in_son[0] = -1.0 + 2.0*s[0];
1565 s_in_son[1] = 1.0 + 2.0*s[1];
1572 s_in_son[0] = -1.0 + 2.0*s[0];
1573 s_in_son[1] = -1.0 + 2.0*s[1];
1578 Node* neighbour_son_node_pt =
1580 get_node_at_local_coordinate(s_in_son);
1583 if(neighbour_son_node_pt!=0)
1587 if(in_neighbouring_tree)
1591 quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1594 return neighbour_son_node_pt;
1615 template<
unsigned INITIAL_NNODE_1D>
1623 if (adopted_father_pt!=0)
1626 father_pt =
dynamic_cast<QuadTree*
>(adopted_father_pt);
1629 else if (Tree_pt!=0)
1632 father_pt =
dynamic_cast<QuadTree*
>(quadtree_pt()->father_pt());
1643 if (father_pt!=0 || initial_p_order!=0)
1650 if (father_el_pt!=0)
1652 unsigned father_p_order = father_el_pt->
p_order();
1654 P_order = father_p_order;
1659 P_order = initial_p_order;
1664 unsigned new_n_node = P_order*P_order;
1667 this->set_n_node(new_n_node);
1670 delete this->integral_pt();
1692 std::ostringstream error_message;
1693 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
1694 error_message <<
"\n integration scheme.\n";
1696 OOMPH_CURRENT_FUNCTION,
1697 OOMPH_EXCEPTION_LOCATION);
1706 template<
unsigned INITIAL_NNODE_1D>
1711 using namespace QuadTreeNames;
1714 unsigned n_p = nnode_1d();
1717 if(Father_bound[n_p].nrow()==0) {setup_father_bounds();}
1723 int son_type = Tree_pt->
son_type();
1732 "Something fishy here: I have no father and yet \n";
1734 "I have no nodes. Who has created me then?!\n";
1737 OOMPH_CURRENT_FUNCTION,
1738 OOMPH_EXCEPTION_LOCATION);
1751 unsigned ntstorage=time_stepper_pt->ntstorage();
1797 for(
unsigned i=0;
i<2;
i++)
1810 if(father_el_pt->
node_pt(0)==0)
1813 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1814 OOMPH_CURRENT_FUNCTION,
1815 OOMPH_EXCEPTION_LOCATION);
1825 for(
unsigned i0=0;i0<n_p;i0++)
1828 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
1830 s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
1832 for(
unsigned i1=0;i1<n_p;i1++)
1835 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
1837 s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
1865 if(created_node_pt!=0)
1868 this->node_pt(jnod) = created_node_pt;
1876 for(
unsigned t=0;
t<ntstorage;
t++)
1885 unsigned n_val_at_node = created_node_pt->
nvalue();
1886 unsigned n_val_from_function = prev_values.size();
1888 unsigned n_var = n_val_at_node < n_val_from_function ?
1889 n_val_at_node : n_val_from_function;
1891 for(
unsigned k=0;k<n_var;k++)
1893 created_node_pt->
set_value(
t,k,prev_values[k]);
1934 template<
unsigned INITIAL_NNODE_1D>
1937 Mesh*
const &mesh_pt,
1954 unsigned ngeom_object=m_el_pt->ngeom_object();
1956 for (
unsigned i=0;
i<ngeom_object;
i++)
1968 "Cloned copy must be a PRefineableQElement<2,INITIAL_NNODE_1D>!",
1969 OOMPH_CURRENT_FUNCTION,
1970 OOMPH_EXCEPTION_LOCATION);
1977 bool clone_is_ok =
true;
1980 clone_is_ok = clone_is_ok && (clone_el_pt->
p_order() == this->p_order());
1984 std::ostringstream err_stream;
1985 err_stream <<
"Clone element has a different p-order from me," 1987 <<
"but it is supposed to be a copy!" << std::endl;
1990 OOMPH_CURRENT_FUNCTION,
1991 OOMPH_EXCEPTION_LOCATION);
1995 clone_is_ok = clone_is_ok && (clone_el_pt->
nnode() == this->nnode());
1999 std::ostringstream err_stream;
2000 err_stream <<
"Clone element has a different number of nodes from me," 2002 <<
"but it is supposed to be a copy!" << std::endl;
2005 OOMPH_CURRENT_FUNCTION,
2006 OOMPH_EXCEPTION_LOCATION);
2010 for(
unsigned n=0; n<this->nnode(); n++)
2012 clone_is_ok = clone_is_ok && (clone_el_pt->
node_pt(n) == this->node_pt(n));
2017 std::ostringstream err_stream;
2018 err_stream <<
"The nodes belonging to the clone element are different" 2020 <<
"from mine, but it is supposed to be a copy!" << std::endl;
2023 OOMPH_CURRENT_FUNCTION,
2024 OOMPH_EXCEPTION_LOCATION);
2040 clone_is_ok = clone_is_ok && (geom_object_pt.size() == clone_geom_object_pt.size());
2041 for(
unsigned n=0; n<std::min(geom_object_pt.size(),clone_geom_object_pt.size()); n++)
2043 clone_is_ok = clone_is_ok && (clone_geom_object_pt[n] == geom_object_pt[n]);
2048 std::ostringstream err_stream;
2049 err_stream <<
"The clone element has different geometric objects" 2051 <<
"from mine, but it is supposed to be a copy!" << std::endl;
2054 "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2055 OOMPH_EXCEPTION_LOCATION);
2069 "Can't handle generalised nodal positions (yet).",
2070 OOMPH_CURRENT_FUNCTION,
2071 OOMPH_EXCEPTION_LOCATION);
2076 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
2079 unsigned ntstorage = time_stepper_pt->
ntstorage();
2085 delete this->integral_pt();
2107 std::ostringstream error_message;
2108 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
2109 error_message <<
"\n integration scheme.\n";
2111 OOMPH_CURRENT_FUNCTION,
2112 OOMPH_EXCEPTION_LOCATION);
2116 this->set_n_node(P_order*P_order);
2137 for(
unsigned i0=0;i0<P_order;i0++)
2140 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
2142 s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
2144 for(
unsigned i1=0;i1<P_order;i1++)
2147 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
2149 s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
2152 jnod= i0 + P_order*i1;
2156 bool node_done=
false;
2164 if (created_node_pt!=0)
2167 this->node_pt(jnod) = created_node_pt;
2175 for(
unsigned t=0;
t<ntstorage;
t++)
2184 unsigned n_val_at_node = created_node_pt->
nvalue();
2185 unsigned n_val_from_function = prev_values.size();
2187 unsigned n_var = n_val_at_node < n_val_from_function ?
2188 n_val_at_node : n_val_from_function;
2190 for(
unsigned k=0;k<n_var;k++)
2192 created_node_pt->
set_value(
t,k,prev_values[k]);
2208 bool is_periodic =
false;
2210 node_created_by_neighbour(s_fraction,is_periodic);
2213 if (created_node_pt!=0)
2222 Node* neighbour_node_pt = created_node_pt;
2236 std::set<unsigned> boundaries;
2246 if (boundaries.size()>1)
2249 "boundaries.size()!=1 seems a bit strange..\n",
2250 OOMPH_CURRENT_FUNCTION,
2251 OOMPH_EXCEPTION_LOCATION);
2255 if(boundaries.size() == 0)
2257 std::ostringstream error_stream;
2259 <<
"Periodic node is not on a boundary...\n" 2261 << created_node_pt->
x(0) <<
" " 2262 << created_node_pt->
x(1) <<
"\n";
2265 OOMPH_CURRENT_FUNCTION,
2266 OOMPH_EXCEPTION_LOCATION);
2272 this->construct_boundary_node(jnod,time_stepper_pt);
2275 make_periodic(neighbour_node_pt);
2278 for (
unsigned t=0;
t<ntstorage;
t++)
2286 clone_el_pt->
get_x(
t,s,x_prev);
2288 for(
unsigned i=0;
i<2;
i++)
2290 created_node_pt->x(
t,
i) = x_prev[
i];
2299 for(std::set<unsigned>::iterator it = boundaries.begin();
2300 it != boundaries.end(); ++it)
2316 created_node_pt->set_coordinates_on_boundary(*it,zeta);
2328 this->node_pt(jnod) = created_node_pt;
2342 bool is_periodic =
false;
2344 node_created_by_son_of_neighbour(s_fraction,is_periodic);
2347 if (created_node_pt!=0)
2356 Node* neighbour_node_pt = created_node_pt;
2370 std::set<unsigned> boundaries;
2380 if (boundaries.size()>1)
2383 "boundaries.size()!=1 seems a bit strange..\n",
2384 OOMPH_CURRENT_FUNCTION,
2385 OOMPH_EXCEPTION_LOCATION);
2389 if(boundaries.size() == 0)
2391 std::ostringstream error_stream;
2393 <<
"Periodic node is not on a boundary...\n" 2395 << created_node_pt->
x(0) <<
" " 2396 << created_node_pt->
x(1) <<
"\n";
2399 OOMPH_CURRENT_FUNCTION,
2400 OOMPH_EXCEPTION_LOCATION);
2406 this->construct_boundary_node(jnod,time_stepper_pt);
2409 make_periodic(neighbour_node_pt);
2412 for (
unsigned t=0;
t<ntstorage;
t++)
2420 clone_el_pt->
get_x(
t,s,x_prev);
2422 for(
unsigned i=0;
i<2;
i++)
2424 created_node_pt->x(
t,
i) = x_prev[
i];
2433 for(std::set<unsigned>::iterator it = boundaries.begin();
2434 it != boundaries.end(); ++it)
2450 created_node_pt->set_coordinates_on_boundary(*it,zeta);
2462 this->node_pt(jnod) = created_node_pt;
2489 std::set<unsigned> boundaries;
2499 if (boundaries.size()>1)
2502 "boundaries.size()!=1 seems a bit strange..\n",
2503 OOMPH_CURRENT_FUNCTION,
2504 OOMPH_EXCEPTION_LOCATION);
2510 if(boundaries.size()> 0)
2513 created_node_pt = this->construct_boundary_node(jnod,time_stepper_pt);
2520 Vector<int> bound_cons(this->ncont_interpolated_values());
2521 clone_el_pt->
get_bcs(my_bound,bound_cons);
2524 unsigned n_value = created_node_pt->
nvalue();
2525 for(
unsigned k=0;k<n_value;k++)
2527 if (bound_cons[k]) {created_node_pt->
pin(k);}
2533 dynamic_cast<SolidNode*
>(created_node_pt);
2534 if (solid_node_pt!=0)
2537 unsigned n_dim = created_node_pt->
ndim();
2542 if (clone_solid_el_pt==0)
2545 "We have a SolidNode outside a refineable SolidElement\n";
2547 "during mesh refinement -- this doesn't make sense";
2551 OOMPH_CURRENT_FUNCTION,
2552 OOMPH_EXCEPTION_LOCATION);
2556 get_solid_bcs(my_bound,solid_bound_cons);
2559 for(
unsigned k=0;k<n_dim;k++)
2561 if (solid_bound_cons[k]) {solid_node_pt->
pin_position(k);}
2572 for(std::set<unsigned>::iterator it = boundaries.begin();
2573 it != boundaries.end(); ++it)
2599 created_node_pt = this->construct_node(jnod,time_stepper_pt);
2615 for (
unsigned t=0;
t<ntstorage;
t++)
2623 clone_el_pt->
get_x(
t,s,x_prev);
2626 for(
unsigned i=0;
i<2;
i++)
2628 created_node_pt->
x(
t,
i) = x_prev[
i];
2633 for (
unsigned t=0;
t<ntstorage;
t++)
2640 unsigned n_value = created_node_pt->
nvalue();
2641 for(
unsigned k=0;k<n_value;k++)
2643 created_node_pt->
set_value(
t,k,prev_values[k]);
2660 "Have not implemented p-refinement for";
2662 "Algebraic p-refineable elements yet\n";
2666 "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2667 OOMPH_EXCEPTION_LOCATION);
2707 if (clone_m_el_pt!=0)
2721 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
2723 "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
2724 error_message +=
"then I should be too....\n";
2727 OOMPH_CURRENT_FUNCTION,
2728 OOMPH_EXCEPTION_LOCATION);
2739 m_el_pt->set_node_update_info(geom_object_pt);
2786 for(
unsigned i0=0;i0<P_order;i0++)
2789 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
2791 s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
2793 for(
unsigned i1=0;i1<P_order;i1++)
2796 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
2798 s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
2801 jnod= i0 + P_order*i1;
2804 for (
unsigned t=0;
t<ntstorage;
t++)
2812 this->get_x(
t,s,x_prev);
2815 for(
unsigned i=0;
i<2;
i++)
2817 this->node_pt(jnod)->x(
t,
i) = x_prev[
i];
2828 this->further_build();
2834 template<
unsigned INITIAL_NNODE_1D>
2848 for(
unsigned i=0;
i<2;
i++)
2850 for(
unsigned j=0;j<2;j++)
2852 psi(2*
i + j) = psi2[
i]*psi1[j];
2865 for(
unsigned i=0;
i<3;
i++)
2867 for(
unsigned j=0;j<3;j++)
2869 psi(3*
i + j) = psi2[
i]*psi1[j];
2882 for(
unsigned i=0;
i<4;
i++)
2884 for(
unsigned j=0;j<4;j++)
2886 psi(4*
i + j) = psi2[
i]*psi1[j];
2899 for(
unsigned i=0;
i<5;
i++)
2901 for(
unsigned j=0;j<5;j++)
2903 psi(5*
i + j) = psi2[
i]*psi1[j];
2916 for(
unsigned i=0;
i<6;
i++)
2918 for(
unsigned j=0;j<6;j++)
2920 psi(6*
i + j) = psi2[
i]*psi1[j];
2933 for(
unsigned i=0;
i<7;
i++)
2935 for(
unsigned j=0;j<7;j++)
2937 psi(7*
i + j) = psi2[
i]*psi1[j];
2943 std::ostringstream error_message;
2944 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
2945 error_message <<
"\n polynomial order for shape functions.\n";
2947 OOMPH_CURRENT_FUNCTION,
2948 OOMPH_EXCEPTION_LOCATION);
2956 template<
unsigned INITIAL_NNODE_1D>
2972 for(
unsigned i=0;
i<2;
i++)
2974 for(
unsigned j=0;j<2;j++)
2977 dpsids(index,0) = psi2[
i]*dpsi1ds[j];
2978 dpsids(index,1) = dpsi2ds[
i]*psi1[j];
2979 psi[index] = psi2[
i]*psi1[j];
2996 for(
unsigned i=0;
i<3;
i++)
2998 for(
unsigned j=0;j<3;j++)
3001 dpsids(index,0) = psi2[
i]*dpsi1ds[j];
3002 dpsids(index,1) = dpsi2ds[
i]*psi1[j];
3003 psi[index] = psi2[
i]*psi1[j];
3020 for(
unsigned i=0;
i<4;
i++)
3022 for(
unsigned j=0;j<4;j++)
3025 dpsids(index,0) = psi2[
i]*dpsi1ds[j];
3026 dpsids(index,1) = dpsi2ds[
i]*psi1[j];
3027 psi[index] = psi2[
i]*psi1[j];
3044 for(
unsigned i=0;
i<5;
i++)
3046 for(
unsigned j=0;j<5;j++)
3049 dpsids(index,0) = psi2[
i]*dpsi1ds[j];
3050 dpsids(index,1) = dpsi2ds[
i]*psi1[j];
3051 psi[index] = psi2[
i]*psi1[j];
3068 for(
unsigned i=0;
i<6;
i++)
3070 for(
unsigned j=0;j<6;j++)
3073 dpsids(index,0) = psi2[
i]*dpsi1ds[j];
3074 dpsids(index,1) = dpsi2ds[
i]*psi1[j];
3075 psi[index] = psi2[
i]*psi1[j];
3092 for(
unsigned i=0;
i<7;
i++)
3094 for(
unsigned j=0;j<7;j++)
3097 dpsids(index,0) = psi2[
i]*dpsi1ds[j];
3098 dpsids(index,1) = dpsi2ds[
i]*psi1[j];
3099 psi[index] = psi2[
i]*psi1[j];
3107 std::ostringstream error_message;
3108 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
3109 error_message <<
"\n polynomial order for shape functions.\n";
3111 OOMPH_CURRENT_FUNCTION,
3112 OOMPH_EXCEPTION_LOCATION);
3121 template<
unsigned INITIAL_NNODE_1D>
3126 std::ostringstream error_message;
3127 error_message <<
"\nd2shape_local currently not implemented for this element\n";
3129 OOMPH_CURRENT_FUNCTION,
3130 OOMPH_EXCEPTION_LOCATION);
3137 template<
unsigned INITIAL_NNODE_1D>
3141 using namespace QuadTreeNames;
3144 unsigned n_sons = this->tree_pt()->nsons();
3146 unsigned max_son_p_order = 0;
3147 for (
unsigned ison=0;ison<n_sons;ison++)
3150 son_p_order[ison] = el_pt->
p_order();
3151 if (son_p_order[ison] > max_son_p_order) max_son_p_order = son_p_order[ison];
3154 unsigned old_Nnode = this->nnode();
3155 unsigned old_P_order = this->p_order();
3157 this->p_order() = max_son_p_order;
3160 delete this->integral_pt();
3161 switch(this->p_order())
3182 std::ostringstream error_message;
3183 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
3184 error_message <<
"\n integration scheme.\n";
3186 OOMPH_CURRENT_FUNCTION,
3187 OOMPH_EXCEPTION_LOCATION);
3192 for (
unsigned n=0; n<old_Nnode; n++)
3194 old_node_pt[n] = this->node_pt(n);
3198 this->set_n_node(this->p_order()*this->p_order());
3201 this->node_pt(0) = old_node_pt[0];
3202 this->node_pt(this->p_order()-1) = old_node_pt[old_P_order-1];
3203 this->node_pt(this->p_order()*(this->p_order()-1))
3204 = old_node_pt[(old_P_order)*(old_P_order-1)];
3205 this->node_pt(this->p_order()*this->p_order()-1)
3206 = old_node_pt[(old_P_order)*(old_P_order)-1];
3209 if(this->p_order() % 2 == 1)
3212 unsigned midpoint = (this->p_order()-1)/2;
3215 this->node_pt(midpoint)
3217 (quadtree_pt()->son_pt(
SW)->object_pt())->vertex_node_pt(1);
3219 this->node_pt(midpoint*this->p_order())
3221 (quadtree_pt()->son_pt(
SW)->object_pt())->vertex_node_pt(2);
3223 this->node_pt((this->p_order()-1)*this->p_order()+midpoint)
3225 (quadtree_pt()->son_pt(
NE)->object_pt())->vertex_node_pt(2);
3227 this->node_pt((midpoint+1)*this->p_order()-1)
3229 (quadtree_pt()->son_pt(
NE)->object_pt())->vertex_node_pt(1);
3237 if(this->node_pt(0)==0)
3240 OOMPH_CURRENT_FUNCTION,
3241 OOMPH_EXCEPTION_LOCATION);
3244 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
3245 unsigned ntstorage = time_stepper_pt->
ntstorage();
3250 unsigned n_p = this->nnode_1d();
3251 for(
unsigned i0=0;i0<n_p;i0++)
3254 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
3256 s[0] = -1.0 + 2.0*s_fraction[0];
3258 for(
unsigned i1=0;i1<n_p;i1++)
3261 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
3263 s[1] = -1.0 + 2.0*s_fraction[1];
3270 bool node_done=
false;
3274 Node* created_node_pt = this->get_node_at_local_coordinate(s);
3278 if (created_node_pt!=0)
3281 this->node_pt(jnod) = created_node_pt;
3295 bool is_periodic =
false;
3297 node_created_by_neighbour(s_fraction,is_periodic);
3300 if(created_node_pt!=0)
3306 "Cannot handle periodic nodes yet",
3307 OOMPH_CURRENT_FUNCTION,
3308 OOMPH_EXCEPTION_LOCATION);
3313 this->node_pt(jnod) = created_node_pt;
3327 using namespace QuadTreeNames;
3330 if(s_fraction[0] < 0.5)
3333 if(s_fraction[1] < 0.5)
3337 s_in_son[0] = -1.0 + 4.0*s_fraction[0];
3338 s_in_son[1] = -1.0 + 4.0*s_fraction[1];
3345 s_in_son[0] = -1.0 + 4.0*s_fraction[0];
3346 s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
3352 if(s_fraction[1] < 0.5)
3356 s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
3357 s_in_son[1] = -1.0 + 4.0*s_fraction[1];
3364 s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
3365 s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
3373 this->tree_pt()->son_pt(son)->object_pt());
3380 if(i0==0) {boundary =
W;}
3382 else if(i0==n_p-1) {boundary =
E;}
3421 std::set<unsigned> boundaries;
3435 if(boundaries.size()>0)
3438 created_node_pt = this->construct_boundary_node(jnod,time_stepper_pt);
3441 Vector<int> bound_cons(this->ncont_interpolated_values());
3442 son_el_pt->
get_bcs(boundary,bound_cons);
3445 unsigned nval = created_node_pt->
nvalue();
3446 for(
unsigned k=0;k<nval;k++)
3448 if(bound_cons[k]) {created_node_pt->
pin(k);}
3454 dynamic_cast<SolidNode*
>(created_node_pt);
3455 if (solid_node_pt!=0)
3458 unsigned n_dim = created_node_pt->
ndim();
3463 if (son_solid_el_pt==0)
3466 "We have a SolidNode outside a refineable SolidElement\n";
3468 "during mesh refinement -- this doesn't make sense\n";
3472 OOMPH_CURRENT_FUNCTION,
3473 OOMPH_EXCEPTION_LOCATION);
3477 get_solid_bcs(boundary,solid_bound_cons);
3480 for(
unsigned k=0;k<n_dim;k++)
3482 if (solid_bound_cons[k]) {solid_node_pt->
pin_position(k);}
3490 for(std::set<unsigned>::iterator it = boundaries.begin();
3491 it != boundaries.end(); ++it)
3515 created_node_pt = this->construct_node(jnod,time_stepper_pt);
3532 for(
unsigned t=0;
t<ntstorage;
t++)
3534 using namespace QuadTreeNames;
3539 son_el_pt->
get_x(
t,s_in_son,x_prev);
3540 for(
unsigned i=0;
i<2;
i++)
3542 created_node_pt->
x(
t,
i) = x_prev[
i];
3548 for(
unsigned t=0;
t<ntstorage;
t++)
3556 for(
unsigned k=0;k<created_node_pt->
nvalue();k++)
3558 created_node_pt->
set_value(
t,k,prev_values[k]);
3573 "Have not implemented rebuilding from sons for";
3575 "Algebraic p-refineable elements yet\n";
3579 "PRefineableQElement<2,INITIAL_NNODE_1D>::rebuild_from_sons()",
3580 OOMPH_EXCEPTION_LOCATION);
3597 for (
unsigned j=0;j<this->nnode();j++)
3601 this->local_coordinate_of_node(j,s_in_node_update_element);
3608 set_node_update_info(
this,s_in_node_update_element,geom_object_pt);
3620 for(
unsigned i0=0;i0<n_p;i0++)
3623 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
3625 s[0] = -1.0 + 2.0*s_fraction[0];
3627 for(
unsigned i1=0;i1<n_p;i1++)
3630 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
3632 s[1] = -1.0 + 2.0*s_fraction[1];
3638 for (
unsigned t=0;
t<ntstorage;
t++)
3646 this->get_x(
t,s,x_prev);
3649 for(
unsigned i=0;
i<2;
i++)
3651 this->node_pt(jnod)->x(
t,
i) = x_prev[
i];
3667 template<
unsigned INITIAL_NNODE_1D>
3681 if(this->macro_elem_pt()!=0)
3687 using namespace QuadTreeNames;
3690 unsigned n_p=nnode_1d();
3699 double max_error_val=0.0;
3702 edges[0] =
S; edges[1] =
N; edges[2] =
W; edges[3] =
E;
3705 for(
unsigned edge_counter=0;edge_counter<4;edge_counter++)
3709 int neigh_edge,diff_level;
3710 bool in_neighbouring_tree;
3716 s_lo_neigh,s_hi_neigh,
3717 neigh_edge,diff_level,
3718 in_neighbouring_tree);
3725 bool is_periodic=
false;
3726 if(in_neighbouring_tree)
3736 bool exclude_this_edge =
false;
3740 exclude_this_edge =
true;
3742 else if(neigh_pt->
nsons() != 0)
3745 exclude_this_edge =
true;
3749 unsigned my_p_order = this->p_order();
3750 unsigned neigh_p_order =
3752 if(my_p_order != neigh_p_order)
3755 exclude_this_edge =
true;
3762 if(dynamic_cast<FiniteElement*>
3763 (neigh_pt->
object_pt())->macro_elem_pt()!=0)
3770 if(!exclude_this_edge)
3774 for(
unsigned i0=0;i0<n_p;i0++)
3777 Node* local_node_pt=0;
3779 switch(edge_counter)
3783 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
3784 s_fraction[1] = 0.0;
3786 local_node_pt = this->node_pt(i0);
3791 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
3792 s_fraction[1] = 1.0;
3794 local_node_pt = this->node_pt(i0 + n_p*(n_p-1));
3799 s_fraction[0] = 0.0;
3800 s_fraction[1] = this->local_one_d_fraction_of_node(i0,1);
3802 local_node_pt = this->node_pt(n_p*i0);
3807 s_fraction[0] = 1.0;
3808 s_fraction[1] = this->local_one_d_fraction_of_node(i0,1);
3810 local_node_pt = this->node_pt(n_p-1 + n_p*i0);
3817 for(
unsigned i=0;
i<2;
i++)
3820 s[
i] = -1.0 + 2.0*s_fraction[
i];
3822 s_in_neighb[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
3823 (s_hi_neigh[
i] - s_lo_neigh[
i]);
3827 for(
unsigned t=0;
t<n_time;
t++)
3834 if(is_periodic==
false)
3836 for(
int i=0;
i<2;
i++)
3839 double err = std::fabs(local_node_pt->
x(
t,
i) - x_in_neighb[
i]);
3845 << local_node_pt->
x(
t,
i)
3846 <<
" " << x_in_neighb[
i]<< std::endl;
3849 << local_node_pt->
x(1) << std::endl;
3854 if (err>max_error_x[
i]) {max_error_x[
i]=err;}
3863 get_interpolated_values(
t,s_in_neighb,values_in_neighb);
3867 this->get_interpolated_values(
t,
s,values);
3873 for(
unsigned ival=0;ival<num_val;ival++)
3875 double err=std::fabs(values[ival] - values_in_neighb[ival]);
3880 << local_node_pt->
x(1) <<
" \n# " 3881 <<
"erru (S)" << err <<
" " << ival <<
" " 3882 << this->get_node_number(local_node_pt) <<
" " 3884 <<
" " << values_in_neighb[ival] << std::endl;
3887 if (err>max_error_val) {max_error_val=err;}
3897 max_error=max_error_x[0];
3898 if (max_error_x[1]>max_error) max_error=max_error_x[1];
3899 if (max_error_val>max_error) max_error=max_error_val;
3903 oomph_info <<
"\n#------------------------------------ \n#Max error " ;
3905 <<
" " << max_error_x[1]
3906 <<
" " << max_error_val << std::endl;
3907 oomph_info <<
"#------------------------------------ \n " << std::endl;
3957 template<
unsigned INITIAL_NNODE_1D>
3960 std::ofstream& output_hangfile)
3962 using namespace QuadTreeNames;
3967 int neigh_edge,diff_level;
3968 bool in_neighbouring_tree;
3972 neigh_pt=this->quadtree_pt()->
3973 gteq_edge_neighbour(my_edge, translate_s, s_lo_neigh,
3974 s_hi_neigh,neigh_edge,diff_level,in_neighbouring_tree);
3981 bool h_type_slave =
false;
3983 bool p_type_slave =
false;
3992 h_type_slave =
true;
3995 else if(neigh_pt->
nsons()==0)
3999 unsigned my_p_order =
4002 unsigned neigh_p_order =
4007 if(neigh_p_order==my_p_order)
4011 else if(neigh_p_order<my_p_order)
4014 p_type_slave =
true;
4036 if (h_type_slave || p_type_slave)
4039 unsigned active_coord_index;
4040 if(my_edge==
N || my_edge==
S) active_coord_index = 0;
4041 else if(my_edge==
E || my_edge==
W) active_coord_index = 1;
4045 "Cannot transform coordinates",
4046 OOMPH_CURRENT_FUNCTION,
4047 OOMPH_EXCEPTION_LOCATION);
4062 const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
4067 bool is_periodic =
false;
4068 if(in_neighbouring_tree)
4069 {is_periodic = this->tree_pt()->root_pt()->is_neighbour_periodic(my_edge);}
4080 "Cannot do mortaring with periodic hanging nodes yet! (Fix me!)",
4081 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4082 OOMPH_EXCEPTION_LOCATION);
4086 unsigned neighbour_node_number=0;
4087 Node* neighbour_node_pt=0;
4090 for(
unsigned i0=0; i0<neigh_n_p; i0++)
4096 neighbour_node_number = i0 + neigh_n_p*(neigh_n_p-1);
4102 neighbour_node_number = i0;
4108 neighbour_node_number = neigh_n_p-1 + neigh_n_p*i0;
4114 neighbour_node_number = neigh_n_p*i0;
4121 OOMPH_CURRENT_FUNCTION,
4122 OOMPH_EXCEPTION_LOCATION);
4126 master_node_number.push_back(neighbour_node_number);
4127 master_node_pt.push_back(neighbour_node_pt);
4131 unsigned local_node_number=0;
4132 Node* local_node_pt=0;
4135 for(
unsigned i0=0; i0<my_n_p; i0++)
4146 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
4147 s_fraction[1] = 1.0;
4148 local_node_number = i0 + my_n_p*(my_n_p-1);
4149 local_node_pt = this->interpolating_node_pt(local_node_number,value_id);
4154 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
4155 s_fraction[1] = 0.0;
4156 local_node_number = i0;
4157 local_node_pt = this->interpolating_node_pt(local_node_number,value_id);
4161 s_fraction[0] = 1.0;
4163 local_one_d_fraction_of_interpolating_node(i0,1,value_id);
4164 local_node_number = my_n_p-1 + my_n_p*i0;
4165 local_node_pt = this->interpolating_node_pt(local_node_number,value_id);
4169 s_fraction[0] = 0.0;
4171 local_one_d_fraction_of_interpolating_node(i0,1,value_id);
4172 local_node_number = my_n_p*i0;
4173 local_node_pt = this->interpolating_node_pt(local_node_number,value_id);
4178 OOMPH_CURRENT_FUNCTION,
4179 OOMPH_EXCEPTION_LOCATION);
4183 slave_node_number.push_back(local_node_number);
4184 slave_node_pt.push_back(local_node_pt);
4187 slave_node_s_fraction.push_back(s_fraction);
4191 const unsigned n_slave_nodes = slave_node_pt.size();
4192 const unsigned n_master_nodes = master_node_pt.size();
4193 const unsigned slave_element_nnode_1d = my_n_p;
4194 const unsigned master_element_nnode_1d = neigh_n_p;
4204 slave_nodal_position,slave_weight);
4208 master_nodal_position,master_weight);
4219 const unsigned n_mortar_vertices = 2;
4222 vertex_pos[1] = this->ninterpolating_node_1d(value_id)-1;
4224 for(
unsigned i=0;
i<my_n_p-n_mortar_vertices;
i++)
4225 {non_vertex_pos[
i] =
i+1;}
4228 for (
unsigned v=0; v<n_mortar_vertices; v++)
4231 bool node_is_shared=
false;
4232 for(
unsigned i=0;
i<master_node_pt.size();
i++)
4234 if(slave_node_pt[vertex_pos[v]]==master_node_pt[
i])
4236 node_is_shared=
true;
4250 for(
unsigned i=0;
i<2;
i++)
4252 s_in_neigh[
i] = s_lo_neigh[
i] +
4253 slave_node_s_fraction[vertex_pos[v]][translate_s[
i]]*
4254 (s_hi_neigh[
i] - s_lo_neigh[
i]);
4262 slave_node_pt[vertex_pos[v]]->set_nonhanging();
4266 for(
unsigned m=0; m<n_master_nodes; m++)
4268 explicit_hang_pt->
set_master_node_pt(m,master_node_pt[m],master_shapes[master_node_number[m]]);
4272 slave_node_pt[vertex_pos[v]]->set_hanging_pt(explicit_hang_pt,-1);
4278 if(n_slave_nodes-n_mortar_vertices > 0)
4286 for (
unsigned i=0;
i<shared_node_M.size();
i++)
4287 {shared_node_M[
i].resize(n_slave_nodes-n_mortar_vertices);}
4290 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
4293 psi[
i] = pow(-1.0,
int((slave_element_nnode_1d-1)-
i-1))
4295 slave_nodal_position[non_vertex_pos[
i]]);
4297 diag_M[
i] = psi[
i]*slave_weight[non_vertex_pos[
i]];
4301 for(
unsigned v=0; v<shared_node_M.size(); v++)
4303 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
4306 if (std::fabs(slave_nodal_position[non_vertex_pos[
i]]
4307 - slave_nodal_position[vertex_pos[v]]) >= 1.0e-8 )
4310 psi[
i] = pow(-1.0,
int((slave_element_nnode_1d-1)-
i-1))
4312 slave_nodal_position[vertex_pos[v]])
4313 / (slave_nodal_position[non_vertex_pos[
i]]
4314 - slave_nodal_position[vertex_pos[v]]);
4318 slave_nodal_position[vertex_pos[v]]))
4322 psi[
i] = pow(-1.0,
int((slave_element_nnode_1d-1)-
i-1))
4324 slave_nodal_position[non_vertex_pos[
i]]);
4331 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4332 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4333 OOMPH_EXCEPTION_LOCATION);
4336 shared_node_M[v][
i] = psi[
i]*slave_weight[vertex_pos[v]];
4345 for(
unsigned i=0;
i<P.size();
i++) {P[
i].resize(n_master_nodes,0.0);}
4359 if (slave_element_nnode_1d >= master_element_nnode_1d)
4362 quadrature_knot = &slave_nodal_position;
4363 quadrature_weight = &slave_weight;
4368 quadrature_knot = &master_nodal_position;
4369 quadrature_weight = &master_weight;
4373 for (
unsigned q=0; q<(*quadrature_weight).size(); q++)
4378 s_on_mortar[0] = (*quadrature_knot)[q];
4381 for(
unsigned k=0; k<n_slave_nodes-n_mortar_vertices; k++)
4384 if (std::fabs(slave_nodal_position[non_vertex_pos[k]]-s_on_mortar[0]) >= 1.0
e-08)
4387 psi[k] = pow(-1.0,
int((slave_element_nnode_1d-1)-k-1))
4389 / (slave_nodal_position[non_vertex_pos[k]]-s_on_mortar[0]);
4396 psi[k] = pow(-1.0,
int((slave_element_nnode_1d-1)-k-1))
4404 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4405 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4406 OOMPH_EXCEPTION_LOCATION);
4412 for(
unsigned i=0;
i<2;
i++)
4414 s_fraction[
i] = (
i==active_coord_index)
4415 ? 0.5*(s_on_mortar[0]+1.0)
4416 : slave_node_s_fraction[vertex_pos[0]][
i];
4421 for(
unsigned i=0;
i<2;
i++)
4423 s_in_neigh[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
4424 (s_hi_neigh[
i] - s_lo_neigh[
i]);
4431 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
4433 for(
unsigned j=0; j<n_master_nodes; j++)
4435 P[
i][j] += master_shapes[master_node_number[j]]*psi[
i]*(*quadrature_weight)[q];
4446 for(
unsigned v=0; v<n_mortar_vertices; v++)
4450 for(
unsigned i=0;
i<2;
i++)
4452 s_fraction[
i] = (
i==active_coord_index)
4453 ? 0.5*(slave_nodal_position[vertex_pos[v]]+1.0)
4454 : slave_node_s_fraction[vertex_pos[0]][
i];
4459 for(
unsigned i=0;
i<2;
i++)
4461 s_in_neigh[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
4462 (s_hi_neigh[
i] - s_lo_neigh[
i]);
4468 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
4470 for(
unsigned k=0; k<n_master_nodes; k++)
4472 P[
i][k] -= master_shapes[master_node_number[k]]*shared_node_M[v][
i];
4479 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
4481 for(
unsigned j=0; j<n_master_nodes; j++)
4490 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
4492 hang_info_pt[
i] =
new HangInfo(n_master_nodes);
4497 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
4499 for(
unsigned j=0; j<n_master_nodes; j++)
4501 hang_info_pt[
i]->set_master_node_pt(j,master_node_pt[j],P[
i][j]);
4507 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
4513 bool node_is_really_shared =
false;
4514 for (
unsigned m=0; m<hang_info_pt[
i]->nmaster(); m++)
4517 if (hang_info_pt[
i]->master_node_pt(m)==slave_node_pt[non_vertex_pos[
i]])
4519 node_is_really_shared =
true;
4523 if(std::fabs(hang_info_pt[
i]->master_weight(m)-1.0) > 1.0e-06)
4527 "Something fishy here -- with shared node at a mortar vertex",
4528 "PRefineableQElemen<2,INITIAL_NNODE_1D>t::quad_hang_helper()",
4529 OOMPH_EXCEPTION_LOCATION);
4536 if(!node_is_really_shared)
4538 slave_node_pt[non_vertex_pos[
i]]->set_hanging_pt(hang_info_pt[
i],-1);
4546 for (
unsigned i=0;
i<n_slave_nodes;
i++)
4549 if (slave_node_pt[
i]->is_hanging())
4552 double weight_sum = 0.0;
4553 for (
unsigned m=0; m<slave_node_pt[
i]->hanging_pt()->nmaster(); m++)
4555 weight_sum += slave_node_pt[
i]->hanging_pt()->master_weight(m);
4559 if (std::fabs(weight_sum-1.0) > 1.0
e-08)
4561 oomph_info <<
"Sum of master weights: " << weight_sum << std::endl;
4563 "Weights in hanging scheme do not sum to 1",
4564 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4565 OOMPH_EXCEPTION_LOCATION);
4572 bool is_master =
false;
4573 for (
unsigned n=0; n<n_master_nodes; n++)
4575 if (slave_node_pt[
i]==master_node_pt[n])
4586 std::ostringstream error_string;
4588 <<
"This node in the slave element is neither" << std::endl
4589 <<
"hanging or shared with a master element." << std::endl;
4593 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4594 OOMPH_EXCEPTION_LOCATION);
4607 for(
unsigned i=0;
i<n_slave_nodes;
i++)
4610 if(slave_node_pt[
i]->is_hanging())
4617 this->local_coordinate_of_node(slave_node_number[
i],s_local);
4622 this->interpolated_x(s_local,x_in_neighb);
4626 slave_node_pt[
i]->x(0)=x_in_neighb[0];
4627 slave_node_pt[
i]->x(1)=x_in_neighb[1];
4640 template<
unsigned INITIAL_NNODE_1D>
4645 unsigned Nnode_1d = this->nnode_1d();
4646 unsigned n0 = n%Nnode_1d;
4647 unsigned n1 = unsigned(
double(n)/
double(Nnode_1d))%Nnode_1d;
4648 unsigned n2 = unsigned(
double(n)/
double(Nnode_1d*Nnode_1d));
4689 std::ostringstream error_message;
4690 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
4691 error_message <<
"\n shape functions.\n";
4693 OOMPH_CURRENT_FUNCTION,
4694 OOMPH_EXCEPTION_LOCATION);
4700 template<
unsigned INITIAL_NNODE_1D>
4704 this->local_coordinate_of_node(n,s_fraction);
4705 s_fraction[0] = 0.5*(s_fraction[0] + 1.0);
4706 s_fraction[1] = 0.5*(s_fraction[1] + 1.0);
4707 s_fraction[2] = 0.5*(s_fraction[2] + 1.0);
4710 template<
unsigned INITIAL_NNODE_1D>
4714 switch(this->nnode_1d())
4735 std::ostringstream error_message;
4736 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
4737 error_message <<
"\n shape functions.\n";
4739 OOMPH_CURRENT_FUNCTION,
4740 OOMPH_EXCEPTION_LOCATION);
4748 template<
unsigned INITIAL_NNODE_1D>
4752 unsigned Nnode_1d = this->nnode_1d();
4759 for (
unsigned i=0;
i<3;
i++)
4764 bool is_found=
false;
4767 if(std::fabs(s[
i] + 1.0) < tol)
4773 else if(std::fabs(s[
i] - 1.0) < tol)
4775 index[
i] = Nnode_1d-1;
4785 for (
unsigned n=1; n<Nnode_1d-1; n++)
4787 if (std::fabs(z[n] - s[
i]) < tol)
4803 return this->node_pt(index[0] + Nnode_1d*index[1] + Nnode_1d*Nnode_1d*index[2]);
4812 template<
unsigned INITIAL_NNODE_1D>
4816 using namespace OcTreeNames;
4822 if(s_fraction[0]==0.0)
4825 if (s_fraction[1]==0.0) {edges.push_back(
LD);}
4826 if (s_fraction[2]==0.0) {edges.push_back(
LB);}
4827 if (s_fraction[1]==1.0) {edges.push_back(
LU);}
4828 if (s_fraction[2]==1.0) {edges.push_back(
LF);}
4831 if(s_fraction[0]==1.0)
4834 if (s_fraction[1]==0.0) {edges.push_back(
RD);}
4835 if (s_fraction[2]==0.0) {edges.push_back(
RB);}
4836 if (s_fraction[1]==1.0) {edges.push_back(
RU);}
4837 if (s_fraction[2]==1.0) {edges.push_back(
RF);}
4840 if(s_fraction[1]==0.0)
4843 if (s_fraction[2]==0.0) {edges.push_back(
DB);}
4844 if (s_fraction[2]==1.0) {edges.push_back(
DF);}
4847 if(s_fraction[1]==1.0)
4850 if (s_fraction[2]==0.0) {edges.push_back(
UB);}
4851 if (s_fraction[2]==1.0) {edges.push_back(
UF);}
4854 if(s_fraction[2]==0.0)
4859 if(s_fraction[2]==1.0)
4865 unsigned n_face = faces.size();
4868 unsigned n_edge = edges.size();
4875 int neigh_face, diff_level;
4879 for(
unsigned j=0;j<n_face;j++)
4883 neigh_pt = octree_pt()->
4884 gteq_face_neighbour(faces[j],translate_s,s_lo_neigh,s_hi_neigh,neigh_face,
4894 bool a_vertex_node_is_built =
false;
4901 "PRefineableQElement<3,INITIAL_NNODE_1D>::node_created_by_neighbour()",
4902 OOMPH_EXCEPTION_LOCATION);
4904 for(
unsigned vnode=0; vnode<neigh_obj_pt->nvertex_node(); vnode++)
4906 if(neigh_obj_pt->vertex_node_pt(vnode)!=0)
4907 a_vertex_node_is_built =
true;
4910 if(a_vertex_node_is_built)
4920 for(
unsigned i=0;
i<3;
i++)
4922 s[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
4923 (s_hi_neigh[
i] - s_lo_neigh[
i]);
4927 Node* neighbour_node_pt =
4931 if(neighbour_node_pt!=0)
4933 return neighbour_node_pt;
4942 for(
unsigned j=0;j<n_edge;j++)
4957 unsigned i_root_edge_neighbour=0;
4960 unsigned nroot_edge_neighbour=0;
4964 bool keep_searching=
true;
4965 while (keep_searching)
4970 neigh_pt = octree_pt()->
4971 gteq_true_edge_neighbour(edges[j],
4972 i_root_edge_neighbour,
4973 nroot_edge_neighbour,
4975 s_lo_neigh,s_hi_neigh,neigh_face,
4985 bool a_vertex_node_is_built =
false;
4992 "PRefineableQElement<3,INITIAL_NNODE_1D>::node_created_by_neighbour()",
4993 OOMPH_EXCEPTION_LOCATION);
4995 for(
unsigned vnode=0; vnode<neigh_obj_pt->nvertex_node(); vnode++)
4997 if(neigh_obj_pt->vertex_node_pt(vnode)!=0)
4998 a_vertex_node_is_built =
true;
5001 if(a_vertex_node_is_built)
5011 for(
unsigned i=0;
i<3;
i++)
5013 s[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
5014 (s_hi_neigh[
i] - s_lo_neigh[
i]);
5018 Node* neighbour_node_pt =
5022 if(neighbour_node_pt!=0)
5024 return neighbour_node_pt;
5031 i_root_edge_neighbour++;
5034 if (i_root_edge_neighbour>=nroot_edge_neighbour)
5036 keep_searching=
false;
5054 template<
unsigned INITIAL_NNODE_1D>
5058 using namespace OcTreeNames;
5064 if(s_fraction[0]==0.0)
5067 if (s_fraction[1]==0.0) {edges.push_back(
LD);}
5068 if (s_fraction[2]==0.0) {edges.push_back(
LB);}
5069 if (s_fraction[1]==1.0) {edges.push_back(
LU);}
5070 if (s_fraction[2]==1.0) {edges.push_back(
LF);}
5073 if(s_fraction[0]==1.0)
5076 if (s_fraction[1]==0.0) {edges.push_back(
RD);}
5077 if (s_fraction[2]==0.0) {edges.push_back(
RB);}
5078 if (s_fraction[1]==1.0) {edges.push_back(
RU);}
5079 if (s_fraction[2]==1.0) {edges.push_back(
RF);}
5082 if(s_fraction[1]==0.0)
5085 if (s_fraction[2]==0.0) {edges.push_back(
DB);}
5086 if (s_fraction[2]==1.0) {edges.push_back(
DF);}
5089 if(s_fraction[1]==1.0)
5092 if (s_fraction[2]==0.0) {edges.push_back(
UB);}
5093 if (s_fraction[2]==1.0) {edges.push_back(
UF);}
5096 if(s_fraction[2]==0.0)
5101 if(s_fraction[2]==1.0)
5107 unsigned n_face = faces.size();
5110 unsigned n_edge = edges.size();
5117 int neigh_face, diff_level;
5121 for(
unsigned j=0;j<n_face;j++)
5125 neigh_pt = octree_pt()->
5126 gteq_face_neighbour(faces[j],translate_s,s_lo_neigh,s_hi_neigh,neigh_face,
5144 for(
unsigned i=0;
i<3;
i++)
5146 s[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
5147 (s_hi_neigh[
i] - s_lo_neigh[
i]);
5151 if(neigh_pt->
nsons()!=0)
5169 s_in_son[0] = 1.0 + 2.0*s[0];
5170 s_in_son[1] = 1.0 + 2.0*s[1];
5171 s_in_son[2] = 1.0 + 2.0*s[2];
5178 s_in_son[0] = 1.0 + 2.0*s[0];
5179 s_in_son[1] = 1.0 + 2.0*s[1];
5180 s_in_son[2] = -1.0 + 2.0*s[2];
5191 s_in_son[0] = 1.0 + 2.0*s[0];
5192 s_in_son[1] = -1.0 + 2.0*s[1];
5193 s_in_son[2] = 1.0 + 2.0*s[2];
5200 s_in_son[0] = 1.0 + 2.0*s[0];
5201 s_in_son[1] = -1.0 + 2.0*s[1];
5202 s_in_son[2] = -1.0 + 2.0*s[2];
5217 s_in_son[0] = -1.0 + 2.0*s[0];
5218 s_in_son[1] = 1.0 + 2.0*s[1];
5219 s_in_son[2] = 1.0 + 2.0*s[2];
5226 s_in_son[0] = -1.0 + 2.0*s[0];
5227 s_in_son[1] = 1.0 + 2.0*s[1];
5228 s_in_son[2] = -1.0 + 2.0*s[2];
5239 s_in_son[0] = -1.0 + 2.0*s[0];
5240 s_in_son[1] = -1.0 + 2.0*s[1];
5241 s_in_son[2] = 1.0 + 2.0*s[2];
5248 s_in_son[0] = -1.0 + 2.0*s[0];
5249 s_in_son[1] = -1.0 + 2.0*s[1];
5250 s_in_son[2] = -1.0 + 2.0*s[2];
5256 Node* neighbour_son_node_pt =
5258 get_node_at_local_coordinate(s_in_son);
5261 if(neighbour_son_node_pt!=0)
5264 return neighbour_son_node_pt;
5274 for(
unsigned j=0;j<n_edge;j++)
5289 unsigned i_root_edge_neighbour=0;
5292 unsigned nroot_edge_neighbour=0;
5296 bool keep_searching=
true;
5297 while (keep_searching)
5302 neigh_pt = octree_pt()->
5303 gteq_true_edge_neighbour(edges[j],
5304 i_root_edge_neighbour,
5305 nroot_edge_neighbour,
5307 s_lo_neigh,s_hi_neigh,neigh_face,
5325 for(
unsigned i=0;
i<3;
i++)
5327 s[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
5328 (s_hi_neigh[
i] - s_lo_neigh[
i]);
5332 if(neigh_pt->
nsons()!=0)
5350 s_in_son[0] = 1.0 + 2.0*s[0];
5351 s_in_son[1] = 1.0 + 2.0*s[1];
5352 s_in_son[2] = 1.0 + 2.0*s[2];
5359 s_in_son[0] = 1.0 + 2.0*s[0];
5360 s_in_son[1] = 1.0 + 2.0*s[1];
5361 s_in_son[2] = -1.0 + 2.0*s[2];
5372 s_in_son[0] = 1.0 + 2.0*s[0];
5373 s_in_son[1] = -1.0 + 2.0*s[1];
5374 s_in_son[2] = 1.0 + 2.0*s[2];
5381 s_in_son[0] = 1.0 + 2.0*s[0];
5382 s_in_son[1] = -1.0 + 2.0*s[1];
5383 s_in_son[2] = -1.0 + 2.0*s[2];
5398 s_in_son[0] = -1.0 + 2.0*s[0];
5399 s_in_son[1] = 1.0 + 2.0*s[1];
5400 s_in_son[2] = 1.0 + 2.0*s[2];
5407 s_in_son[0] = -1.0 + 2.0*s[0];
5408 s_in_son[1] = 1.0 + 2.0*s[1];
5409 s_in_son[2] = -1.0 + 2.0*s[2];
5420 s_in_son[0] = -1.0 + 2.0*s[0];
5421 s_in_son[1] = -1.0 + 2.0*s[1];
5422 s_in_son[2] = 1.0 + 2.0*s[2];
5429 s_in_son[0] = -1.0 + 2.0*s[0];
5430 s_in_son[1] = -1.0 + 2.0*s[1];
5431 s_in_son[2] = -1.0 + 2.0*s[2];
5437 Node* neighbour_son_node_pt =
5439 get_node_at_local_coordinate(s_in_son);
5442 if(neighbour_son_node_pt!=0)
5445 return neighbour_son_node_pt;
5453 i_root_edge_neighbour++;
5456 if (i_root_edge_neighbour>=nroot_edge_neighbour)
5458 keep_searching=
false;
5475 template<
unsigned INITIAL_NNODE_1D>
5483 if (adopted_father_pt!=0)
5486 father_pt =
dynamic_cast<OcTree*
>(adopted_father_pt);
5489 else if (Tree_pt!=0)
5492 father_pt =
dynamic_cast<OcTree*
>(octree_pt()->father_pt());
5503 if (father_pt!=0 || initial_p_order!=0)
5510 if (father_el_pt!=0)
5512 unsigned father_p_order = father_el_pt->
p_order();
5514 P_order = father_p_order;
5519 P_order = initial_p_order;
5524 unsigned new_n_node = P_order*P_order*P_order;
5527 this->set_n_node(new_n_node);
5530 delete this->integral_pt();
5552 std::ostringstream error_message;
5553 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
5554 error_message <<
"\n integration scheme.\n";
5556 OOMPH_CURRENT_FUNCTION,
5557 OOMPH_EXCEPTION_LOCATION);
5566 template<
unsigned INITIAL_NNODE_1D>
5571 using namespace OcTreeNames;
5574 unsigned n_p = nnode_1d();
5577 if(Father_bound[n_p].nrow()==0) {setup_father_bounds();}
5580 OcTree* father_pt =
dynamic_cast<OcTree*
>(octree_pt()->father_pt());
5583 int son_type = Tree_pt->
son_type();
5592 "Something fishy here: I have no father and yet \n";
5594 "I have no nodes. Who has created me then?!\n";
5597 "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
5598 OOMPH_EXCEPTION_LOCATION);
5611 unsigned ntstorage=time_stepper_pt->ntstorage();
5717 if(father_el_pt->
node_pt(0)==0)
5720 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
5721 "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
5722 OOMPH_EXCEPTION_LOCATION);
5730 for(
unsigned i0=0;i0<n_p;i0++)
5733 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
5735 s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
5737 for(
unsigned i1=0;i1<n_p;i1++)
5740 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
5742 s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
5744 for(
unsigned i2=0;i2<n_p;i2++)
5747 s_fraction[2] = this->local_one_d_fraction_of_node(i2,2);
5749 s[2] = s_lo[2] + (s_hi[2]-s_lo[2])*s_fraction[2];
5752 jnod= i0 + n_p*i1 + n_p*n_p*i2;
5760 if(created_node_pt!=0)
5763 this->node_pt(jnod) = created_node_pt;
5771 for(
unsigned t=0;
t<ntstorage;
t++)
5780 unsigned n_val_at_node = created_node_pt->
nvalue();
5781 unsigned n_val_from_function = prev_values.size();
5783 unsigned n_var = n_val_at_node < n_val_from_function ?
5784 n_val_at_node : n_val_from_function;
5786 for(
unsigned k=0;k<n_var;k++)
5788 created_node_pt->
set_value(
t,k,prev_values[k]);
5805 template<
unsigned INITIAL_NNODE_1D>
5808 Mesh*
const &mesh_pt,
5819 "Cloned copy must be a PRefineableQElement<3,INITIAL_NNODE_1D>!",
5820 OOMPH_CURRENT_FUNCTION,
5821 OOMPH_EXCEPTION_LOCATION);
5828 bool clone_is_ok =
true;
5831 clone_is_ok = clone_is_ok && (clone_el_pt->
p_order() == this->p_order());
5835 std::ostringstream err_stream;
5836 err_stream <<
"Clone element has a different p-order from me," 5838 <<
"but it is supposed to be a copy!" << std::endl;
5841 OOMPH_CURRENT_FUNCTION,
5842 OOMPH_EXCEPTION_LOCATION);
5846 clone_is_ok = clone_is_ok && (clone_el_pt->
nnode() == this->nnode());
5850 std::ostringstream err_stream;
5851 err_stream <<
"Clone element has a different number of nodes from me," 5853 <<
"but it is supposed to be a copy!" << std::endl;
5856 OOMPH_CURRENT_FUNCTION,
5857 OOMPH_EXCEPTION_LOCATION);
5861 for(
unsigned n=0; n<this->nnode(); n++)
5863 clone_is_ok = clone_is_ok && (clone_el_pt->
node_pt(n) == this->node_pt(n));
5868 std::ostringstream err_stream;
5869 err_stream <<
"The nodes belonging to the clone element are different" 5871 <<
"from mine, but it is supposed to be a copy!" << std::endl;
5874 OOMPH_CURRENT_FUNCTION,
5875 OOMPH_EXCEPTION_LOCATION);
5888 "Can't handle generalised nodal positions (yet).",
5889 OOMPH_CURRENT_FUNCTION,
5890 OOMPH_EXCEPTION_LOCATION);
5895 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
5898 unsigned ntstorage = time_stepper_pt->
ntstorage();
5901 this->P_order += inc;
5904 delete this->integral_pt();
5905 switch(this->P_order)
5926 std::ostringstream error_message;
5927 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
5928 error_message <<
"\n integration scheme.\n";
5930 OOMPH_CURRENT_FUNCTION,
5931 OOMPH_EXCEPTION_LOCATION);
5935 this->set_n_node(this->P_order*this->P_order*this->P_order);
5958 for(
unsigned i0=0;i0<this->P_order;i0++)
5961 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
5963 s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
5965 for(
unsigned i1=0;i1<this->P_order;i1++)
5968 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
5970 s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
5972 for(
unsigned i2=0;i2<this->P_order;i2++)
5975 s_fraction[2] = this->local_one_d_fraction_of_node(i2,2);
5977 s[2] = s_lo[2] + (s_hi[2]-s_lo[2])*s_fraction[2];
5980 jnod= i0 + this->P_order*i1 + this->P_order*this->P_order*i2;
5984 bool node_done=
false;
5993 if (created_node_pt!=0)
5996 this->node_pt(jnod) = created_node_pt;
6004 for(
unsigned t=0;
t<ntstorage;
t++)
6013 unsigned n_val_at_node = created_node_pt->
nvalue();
6014 unsigned n_val_from_function = prev_values.size();
6016 unsigned n_var = n_val_at_node < n_val_from_function ?
6017 n_val_at_node : n_val_from_function;
6019 for(
unsigned k=0;k<n_var;k++)
6021 created_node_pt->
set_value(
t,k,prev_values[k]);
6038 this->node_created_by_neighbour(s_fraction);
6041 if (created_node_pt!=0)
6043 this->node_pt(jnod) = created_node_pt;
6057 this->node_created_by_son_of_neighbour(s_fraction);
6060 if (created_node_pt!=0)
6062 this->node_pt(jnod) = created_node_pt;
6102 else if(i1==this->P_order-1)
6156 else if(i2==this->P_order-1)
6195 std::set<unsigned> boundaries;
6205 if (boundaries.size()>2)
6208 "boundaries.size()>2 seems a bit strange..\n",
6209 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6210 OOMPH_EXCEPTION_LOCATION);
6216 if(boundaries.size()> 0)
6219 created_node_pt = this->construct_boundary_node(jnod,time_stepper_pt);
6226 Vector<int> bound_cons(this->ncont_interpolated_values());
6227 clone_el_pt->
get_bcs(my_bound,bound_cons);
6230 unsigned n_value = created_node_pt->
nvalue();
6231 for(
unsigned k=0;k<n_value;k++)
6233 if (bound_cons[k]) {created_node_pt->
pin(k);}
6239 dynamic_cast<SolidNode*
>(created_node_pt);
6240 if (solid_node_pt!=0)
6243 unsigned n_dim = created_node_pt->
ndim();
6248 if (clone_solid_el_pt==0)
6251 "We have a SolidNode outside a refineable SolidElement\n";
6253 "during mesh refinement -- this doesn't make sense";
6257 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6258 OOMPH_EXCEPTION_LOCATION);
6262 get_solid_bcs(my_bound,solid_bound_cons);
6265 for(
unsigned k=0;k<n_dim;k++)
6267 if (solid_bound_cons[k]) {solid_node_pt->
pin_position(k);}
6277 for(std::set<unsigned>::iterator it = boundaries.begin();
6278 it != boundaries.end(); ++it)
6304 created_node_pt = this->construct_node(jnod,time_stepper_pt);
6320 for (
unsigned t=0;
t<ntstorage;
t++)
6328 clone_el_pt->
get_x(
t,s,x_prev);
6331 for(
unsigned i=0;
i<3;
i++)
6333 created_node_pt->
x(
t,
i) = x_prev[
i];
6338 for (
unsigned t=0;
t<ntstorage;
t++)
6346 unsigned n_value = created_node_pt->
nvalue();
6347 for(
unsigned k=0;k<n_value;k++)
6349 created_node_pt->
set_value(
t,k,prev_values[k]);
6366 "Have not implemented p-refinement for";
6368 "Algebraic p-refineable elements yet\n";
6372 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6373 OOMPH_EXCEPTION_LOCATION);
6411 for(
unsigned i0=0;i0<this->P_order;i0++)
6414 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
6416 s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
6418 for(
unsigned i1=0;i1<this->P_order;i1++)
6421 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
6423 s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
6425 for(
unsigned i2=0;i2<this->P_order;i2++)
6428 s_fraction[2] = this->local_one_d_fraction_of_node(i2,2);
6430 s[2] = s_lo[2] + (s_hi[2]-s_lo[2])*s_fraction[2];
6433 jnod= i0 + this->P_order*i1 + this->P_order*this->P_order*i2;
6436 for (
unsigned t=0;
t<ntstorage;
t++)
6444 this->get_x(
t,s,x_prev);
6447 for(
unsigned i=0;
i<3;
i++)
6449 this->node_pt(jnod)->x(
t,
i) = x_prev[
i];
6465 if (clone_m_el_pt!=0)
6479 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
6481 "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
6482 error_message +=
"then I should be too....\n";
6485 "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6486 OOMPH_EXCEPTION_LOCATION);
6497 m_el_pt->set_node_update_info(geom_object_pt);
6504 this->further_build();
6510 template<
unsigned INITIAL_NNODE_1D>
6524 for(
unsigned i=0;
i<P_order;
i++)
6526 for(
unsigned j=0;j<P_order;j++)
6528 for(
unsigned k=0;k<P_order;k++)
6530 psi(P_order*P_order*
i + P_order*j + k) = psi3[
i]*psi2[j]*psi1[k];
6544 for(
unsigned i=0;
i<P_order;
i++)
6546 for(
unsigned j=0;j<P_order;j++)
6548 for(
unsigned k=0;k<P_order;k++)
6550 psi(P_order*P_order*
i + P_order*j + k) = psi3[
i]*psi2[j]*psi1[k];
6564 for(
unsigned i=0;
i<P_order;
i++)
6566 for(
unsigned j=0;j<P_order;j++)
6568 for(
unsigned k=0;k<P_order;k++)
6570 psi(P_order*P_order*
i + P_order*j + k) = psi3[
i]*psi2[j]*psi1[k];
6584 for(
unsigned i=0;
i<P_order;
i++)
6586 for(
unsigned j=0;j<P_order;j++)
6588 for(
unsigned k=0;k<P_order;k++)
6590 psi(P_order*P_order*
i + P_order*j + k) = psi3[
i]*psi2[j]*psi1[k];
6604 for(
unsigned i=0;
i<P_order;
i++)
6606 for(
unsigned j=0;j<P_order;j++)
6608 for(
unsigned k=0;k<P_order;k++)
6610 psi(P_order*P_order*
i + P_order*j + k) = psi3[
i]*psi2[j]*psi1[k];
6624 for(
unsigned i=0;
i<P_order;
i++)
6626 for(
unsigned j=0;j<P_order;j++)
6628 for(
unsigned k=0;k<P_order;k++)
6630 psi(P_order*P_order*
i + P_order*j + k) = psi3[
i]*psi2[j]*psi1[k];
6637 std::ostringstream error_message;
6638 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
6639 error_message <<
"\n polynomial order for shape functions.\n";
6641 OOMPH_CURRENT_FUNCTION,
6642 OOMPH_EXCEPTION_LOCATION);
6650 template<
unsigned INITIAL_NNODE_1D>
6666 for(
unsigned i=0;
i<P_order;
i++)
6668 for(
unsigned j=0;j<P_order;j++)
6670 for(
unsigned k=0;k<P_order;k++)
6673 dpsids(index,0) = psi3[
i]*psi2[j]*dpsi1ds[k];
6674 dpsids(index,1) = psi3[
i]*dpsi2ds[j]*psi1[k];
6675 dpsids(index,2) = dpsi3ds[
i]*psi2[j]*psi1[k];
6676 psi[index] = psi3[
i]*psi2[j]*psi1[k];
6694 for(
unsigned i=0;
i<P_order;
i++)
6696 for(
unsigned j=0;j<P_order;j++)
6698 for(
unsigned k=0;k<P_order;k++)
6701 dpsids(index,0) = psi3[
i]*psi2[j]*dpsi1ds[k];
6702 dpsids(index,1) = psi3[
i]*dpsi2ds[j]*psi1[k];
6703 dpsids(index,2) = dpsi3ds[
i]*psi2[j]*psi1[k];
6704 psi[index] = psi3[
i]*psi2[j]*psi1[k];
6722 for(
unsigned i=0;
i<P_order;
i++)
6724 for(
unsigned j=0;j<P_order;j++)
6726 for(
unsigned k=0;k<P_order;k++)
6729 dpsids(index,0) = psi3[
i]*psi2[j]*dpsi1ds[k];
6730 dpsids(index,1) = psi3[
i]*dpsi2ds[j]*psi1[k];
6731 dpsids(index,2) = dpsi3ds[
i]*psi2[j]*psi1[k];
6732 psi[index] = psi3[
i]*psi2[j]*psi1[k];
6750 for(
unsigned i=0;
i<P_order;
i++)
6752 for(
unsigned j=0;j<P_order;j++)
6754 for(
unsigned k=0;k<P_order;k++)
6757 dpsids(index,0) = psi3[
i]*psi2[j]*dpsi1ds[k];
6758 dpsids(index,1) = psi3[
i]*dpsi2ds[j]*psi1[k];
6759 dpsids(index,2) = dpsi3ds[
i]*psi2[j]*psi1[k];
6760 psi[index] = psi3[
i]*psi2[j]*psi1[k];
6778 for(
unsigned i=0;
i<P_order;
i++)
6780 for(
unsigned j=0;j<P_order;j++)
6782 for(
unsigned k=0;k<P_order;k++)
6785 dpsids(index,0) = psi3[
i]*psi2[j]*dpsi1ds[k];
6786 dpsids(index,1) = psi3[
i]*dpsi2ds[j]*psi1[k];
6787 dpsids(index,2) = dpsi3ds[
i]*psi2[j]*psi1[k];
6788 psi[index] = psi3[
i]*psi2[j]*psi1[k];
6806 for(
unsigned i=0;
i<P_order;
i++)
6808 for(
unsigned j=0;j<P_order;j++)
6810 for(
unsigned k=0;k<P_order;k++)
6813 dpsids(index,0) = psi3[
i]*psi2[j]*dpsi1ds[k];
6814 dpsids(index,1) = psi3[
i]*dpsi2ds[j]*psi1[k];
6815 dpsids(index,2) = dpsi3ds[
i]*psi2[j]*psi1[k];
6816 psi[index] = psi3[
i]*psi2[j]*psi1[k];
6825 std::ostringstream error_message;
6826 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
6827 error_message <<
"\n polynomial order for shape functions.\n";
6829 OOMPH_CURRENT_FUNCTION,
6830 OOMPH_EXCEPTION_LOCATION);
6839 template<
unsigned INITIAL_NNODE_1D>
6844 std::ostringstream error_message;
6845 error_message <<
"\nd2shape_local currently not implemented for this element\n";
6847 OOMPH_CURRENT_FUNCTION,
6848 OOMPH_EXCEPTION_LOCATION);
6855 template<
unsigned INITIAL_NNODE_1D>
6860 unsigned n_sons = this->tree_pt()->nsons();
6862 unsigned max_son_p_order = 0;
6863 for (
unsigned ison=0;ison<n_sons;ison++)
6866 son_p_order[ison] = el_pt->
p_order();
6867 if (son_p_order[ison] > max_son_p_order) max_son_p_order = son_p_order[ison];
6870 unsigned old_Nnode = this->nnode();
6871 unsigned old_P_order = this->p_order();
6873 this->p_order() = max_son_p_order;
6876 delete this->integral_pt();
6877 switch(this->p_order())
6898 std::ostringstream error_message;
6899 error_message <<
"\nERROR: Exceeded maximum polynomial order for";
6900 error_message <<
"\n integration scheme.\n";
6902 "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
6903 OOMPH_EXCEPTION_LOCATION);
6908 for (
unsigned n=0; n<old_Nnode; n++)
6910 old_node_pt[n] = this->node_pt(n);
6914 this->set_n_node(this->p_order()*this->p_order()*this->p_order());
6919 this->node_pt(this->p_order()-1)
6920 = old_node_pt[old_P_order-1];
6921 this->node_pt(this->p_order()*(this->p_order()-1))
6922 = old_node_pt[(old_P_order)*(old_P_order-1)];
6923 this->node_pt(this->p_order()*this->p_order()-1)
6924 = old_node_pt[(old_P_order)*(old_P_order)-1];
6925 this->node_pt(this->p_order()*this->p_order()*(this->p_order()-1))
6926 = old_node_pt[old_P_order*old_P_order*(old_P_order-1)];
6927 this->node_pt((this->p_order()*this->p_order()+1)*(this->p_order()-1))
6928 = old_node_pt[(old_P_order*old_P_order+1)*(old_P_order-1)];
6929 this->node_pt(this->p_order()*(this->p_order()+1)*(this->p_order()-1))
6930 = old_node_pt[old_P_order*(old_P_order+1)*(old_P_order-1)];
6931 this->node_pt(this->p_order()*this->p_order()*this->p_order()-1)
6932 = old_node_pt[old_P_order*old_P_order*old_P_order-1];
6935 if(this->p_order() % 2 == 1)
6938 unsigned n_p = this->p_order();
6939 unsigned m = (n_p-1)/2;
6941 unsigned s0space = m;
6942 unsigned s1space = m*n_p;
6943 unsigned s2space = m*n_p*n_p;
6947 this->node_pt(1*s0space + 0*s1space + 0*s2space)
6951 this->node_pt(0*s0space + 1*s1space + 0*s2space)
6955 this->node_pt(1*s0space + 1*s1space + 0*s2space)
6959 this->node_pt(2*s0space + 1*s1space + 0*s2space)
6963 this->node_pt(1*s0space + 2*s1space + 0*s2space)
6969 this->node_pt(0*s0space + 0*s1space + 1*s2space)
6973 this->node_pt(1*s0space + 0*s1space + 1*s2space)
6977 this->node_pt(2*s0space + 0*s1space + 1*s2space)
6981 this->node_pt(0*s0space + 1*s1space + 1*s2space)
6985 this->node_pt(1*s0space + 1*s1space + 1*s2space)
6989 this->node_pt(2*s0space + 1*s1space + 1*s2space)
6993 this->node_pt(0*s0space + 2*s1space + 1*s2space)
6997 this->node_pt(1*s0space + 2*s1space + 1*s2space)
7001 this->node_pt(2*s0space + 2*s1space + 1*s2space)
7007 this->node_pt(1*s0space + 0*s1space + 2*s2space)
7011 this->node_pt(0*s0space + 1*s1space + 2*s2space)
7015 this->node_pt(1*s0space + 1*s1space + 2*s2space)
7019 this->node_pt(2*s0space + 1*s1space + 2*s2space)
7023 this->node_pt(1*s0space + 2*s1space + 2*s2space)
7030 if(this->node_pt(0)==0)
7033 "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
7034 OOMPH_EXCEPTION_LOCATION);
7037 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
7038 unsigned ntstorage = time_stepper_pt->
ntstorage();
7043 unsigned n_p = this->nnode_1d();
7044 for(
unsigned i0=0;i0<n_p;i0++)
7047 s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
7049 s[0] = -1.0 + 2.0*s_fraction[0];
7051 for(
unsigned i1=0;i1<n_p;i1++)
7054 s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
7056 s[1] = -1.0 + 2.0*s_fraction[1];
7058 for(
unsigned i2=0;i2<n_p;i2++)
7061 s_fraction[2] = this->local_one_d_fraction_of_node(i2,2);
7063 s[2] = -1.0 + 2.0*s_fraction[2];
7066 jnod = i0 + n_p*i1 + n_p*n_p*i2;
7070 bool node_done=
false;
7074 Node* created_node_pt = this->get_node_at_local_coordinate(s);
7078 if (created_node_pt!=0)
7081 this->node_pt(jnod) = created_node_pt;
7096 node_created_by_neighbour(s_fraction);
7099 if(created_node_pt!=0)
7101 this->node_pt(jnod) = created_node_pt;
7117 if(s_fraction[0] < 0.5)
7120 if(s_fraction[1] < 0.5)
7123 if(s_fraction[2] < 0.5)
7127 s_in_son[0] = -1.0 + 4.0*s_fraction[0];
7128 s_in_son[1] = -1.0 + 4.0*s_fraction[1];
7129 s_in_son[2] = -1.0 + 4.0*s_fraction[2];
7136 s_in_son[0] = -1.0 + 4.0*s_fraction[0];
7137 s_in_son[1] = -1.0 + 4.0*s_fraction[1];
7138 s_in_son[2] = -1.0 + 4.0*(s_fraction[2]-0.5);
7149 s_in_son[0] = -1.0 + 4.0*s_fraction[0];
7150 s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
7151 s_in_son[2] = -1.0 + 4.0*s_fraction[2];
7158 s_in_son[0] = -1.0 + 4.0*s_fraction[0];
7159 s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
7160 s_in_son[2] = -1.0 + 4.0*(s_fraction[2]-0.5);
7175 s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
7176 s_in_son[1] = -1.0 + 4.0*s_fraction[1];
7177 s_in_son[2] = -1.0 + 4.0*s_fraction[2];
7184 s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
7185 s_in_son[1] = -1.0 + 4.0*s_fraction[1];
7186 s_in_son[2] = -1.0 + 4.0*(s_fraction[2]-0.5);
7197 s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
7198 s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
7199 s_in_son[2] = -1.0 + 4.0*s_fraction[2];
7206 s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
7207 s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
7208 s_in_son[2] = -1.0 + 4.0*(s_fraction[2]-0.5);
7217 this->tree_pt()->son_pt(son)->object_pt());
7337 std::set<unsigned> boundaries;
7349 if(boundaries.size()>0)
7352 created_node_pt = this->construct_boundary_node(jnod,time_stepper_pt);
7355 Vector<int> bound_cons(this->ncont_interpolated_values());
7356 son_el_pt->
get_bcs(boundary,bound_cons);
7359 unsigned nval = created_node_pt->
nvalue();
7360 for(
unsigned k=0;k<nval;k++)
7362 if(bound_cons[k]) {created_node_pt->
pin(k);}
7368 dynamic_cast<SolidNode*
>(created_node_pt);
7369 if (solid_node_pt!=0)
7372 unsigned n_dim = created_node_pt->
ndim();
7377 if (son_solid_el_pt==0)
7380 "We have a SolidNode outside a refineable SolidElement\n";
7382 "during mesh refinement -- this doesn't make sense\n";
7386 "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
7387 OOMPH_EXCEPTION_LOCATION);
7391 get_solid_bcs(boundary,solid_bound_cons);
7394 for(
unsigned k=0;k<n_dim;k++)
7396 if (solid_bound_cons[k]) {solid_node_pt->
pin_position(k);}
7404 for(std::set<unsigned>::iterator it = boundaries.begin();
7405 it != boundaries.end(); ++it)
7429 created_node_pt = this->construct_node(jnod,time_stepper_pt);
7446 for(
unsigned t=0;
t<ntstorage;
t++)
7448 using namespace QuadTreeNames;
7453 son_el_pt->
get_x(
t,s_in_son,x_prev);
7454 for(
unsigned i=0;
i<3;
i++)
7456 created_node_pt->
x(
t,
i) = x_prev[
i];
7462 for(
unsigned t=0;
t<ntstorage;
t++)
7470 for(
unsigned k=0;k<created_node_pt->
nvalue();k++)
7472 created_node_pt->
set_value(
t,k,prev_values[k]);
7487 "Have not implemented rebuilding from sons for";
7489 "Algebraic p-refineable elements yet\n";
7493 "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
7494 OOMPH_EXCEPTION_LOCATION);
7511 template<
unsigned INITIAL_NNODE_1D>
7526 template<
unsigned INITIAL_NNODE_1D>
7529 const int &my_face, std::ofstream& output_hangfile)
7531 using namespace OcTreeNames;
7536 int neigh_face,diff_level;
7540 neigh_pt=this->octree_pt()->
7541 gteq_face_neighbour(my_face, translate_s, s_lo_neigh,
7542 s_hi_neigh,neigh_face,diff_level);
7549 bool h_type_slave =
false;
7551 bool p_type_slave =
false;
7560 h_type_slave =
true;
7563 else if(neigh_pt->
nsons()==0)
7567 unsigned my_p_order =
7570 unsigned neigh_p_order =
7575 if(neigh_p_order==my_p_order)
7579 else if(neigh_p_order<my_p_order)
7582 p_type_slave =
true;
7604 if (h_type_slave || p_type_slave ||
true)
7612 face_edge_other_face[0] =
B;
7614 face_edge_other_face[1] =
F;
7616 face_edge_other_face[2] =
D;
7618 face_edge_other_face[3] =
U;
7622 face_edge_other_face[0] =
B;
7624 face_edge_other_face[1] =
F;
7626 face_edge_other_face[2] =
D;
7628 face_edge_other_face[3] =
U;
7632 face_edge_other_face[0] =
B;
7634 face_edge_other_face[1] =
F;
7636 face_edge_other_face[2] =
L;
7638 face_edge_other_face[3] =
R;
7642 face_edge_other_face[0] =
B;
7644 face_edge_other_face[1] =
F;
7646 face_edge_other_face[2] =
L;
7648 face_edge_other_face[3] =
R;
7652 face_edge_other_face[0] =
D;
7654 face_edge_other_face[1] =
U;
7656 face_edge_other_face[2] =
L;
7658 face_edge_other_face[3] =
R;
7662 face_edge_other_face[0] =
D;
7664 face_edge_other_face[1] =
U;
7666 face_edge_other_face[2] =
L;
7668 face_edge_other_face[3] =
R;
7672 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
7673 OOMPH_EXCEPTION_LOCATION);
7677 for(
unsigned i=0;
i<4;
i++)
7680 unsigned my_edge = face_edge[
i];
7687 int neigh_edge=0,edge_diff_level=0;
7688 unsigned edge_p_order=
7695 int tmp_neigh_edge,tmp_edge_diff_level;
7698 unsigned i_root_edge_neighbour=0;
7701 unsigned nroot_edge_neighbour=0;
7704 bool my_edge_is_master =
true;
7709 bool keep_searching =
true;
7710 bool search_faces =
false;
7711 bool first_face_searched =
false;
7713 while (keep_searching)
7716 OcTree* tmp_edge_neigh_pt;
7722 if(!first_face_searched)
7725 tmp_edge_neigh_pt=this->octree_pt()->
7726 gteq_face_neighbour(my_face,
7727 tmp_edge_translate_s,
7728 tmp_edge_s_lo_neigh, tmp_edge_s_hi_neigh,
7729 tmp_neigh_edge, tmp_edge_diff_level);
7732 first_face_searched =
true;
7737 tmp_edge_neigh_pt=this->octree_pt()->
7738 gteq_face_neighbour(face_edge_other_face[
i],
7739 tmp_edge_translate_s,
7740 tmp_edge_s_lo_neigh, tmp_edge_s_hi_neigh,
7741 tmp_neigh_edge, tmp_edge_diff_level);
7744 keep_searching =
false;
7752 tmp_edge_neigh_pt=this->octree_pt()->
7753 gteq_true_edge_neighbour(my_edge,
7754 i_root_edge_neighbour,
7755 nroot_edge_neighbour,
7756 tmp_edge_translate_s,
7757 tmp_edge_s_lo_neigh, tmp_edge_s_hi_neigh,
7758 tmp_neigh_edge, tmp_edge_diff_level);
7768 bool new_edge_master =
false;
7771 if(tmp_edge_neigh_pt!=0)
7774 if(tmp_edge_diff_level<edge_diff_level)
7778 new_edge_master =
true;
7780 edge_diff_level = tmp_edge_diff_level;
7783 (tmp_edge_neigh_pt->
object_pt())->p_order();
7786 else if(tmp_edge_diff_level==edge_diff_level && tmp_edge_neigh_pt->
nsons()==0)
7793 unsigned tmp_edge_neigh_p_order =
7795 (tmp_edge_neigh_pt->
object_pt())->p_order();
7798 if(tmp_edge_neigh_p_order==edge_p_order)
7802 else if(tmp_edge_neigh_p_order<edge_p_order)
7806 new_edge_master =
true;
7808 edge_diff_level = tmp_edge_diff_level;
7811 (tmp_edge_neigh_pt->
object_pt())->p_order();
7835 edge_neigh_pt = tmp_edge_neigh_pt;
7836 edge_translate_s = tmp_edge_translate_s;
7837 edge_s_lo_neigh = tmp_edge_s_lo_neigh;
7838 edge_s_hi_neigh = tmp_edge_s_hi_neigh;
7839 neigh_edge = tmp_neigh_edge;
7840 edge_diff_level = tmp_edge_diff_level;
7843 my_edge_is_master =
false;
7848 i_root_edge_neighbour++;
7855 if (i_root_edge_neighbour >= nroot_edge_neighbour)
7859 search_faces =
true;
7865 if (!my_edge_is_master)
7868 unsigned active_coord_index;
7875 active_coord_index = 0;
7881 active_coord_index = 1;
7887 active_coord_index = 2;
7891 "Cannot transform coordinates",
7892 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
7893 OOMPH_EXCEPTION_LOCATION);
7908 const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
7912 unsigned neighbour_node_number=0;
7913 Node* neighbour_node_pt=0;
7916 bool master_is_not_edge =
false;
7917 for(
unsigned i0=0; i0<neigh_n_p; i0++)
7919 const unsigned s0space = 1;
7920 const unsigned s1space = neigh_n_p;
7921 const unsigned s2space = neigh_n_p*neigh_n_p;
7927 neighbour_node_number = i0*s0space;
7932 neighbour_node_number = (neigh_n_p-1)*s1space + i0*s0space;
7937 neighbour_node_number = (neigh_n_p-1)*s2space + i0*s0space;
7942 neighbour_node_number = (neigh_n_p-1)*s1space + (neigh_n_p-1)*s2space + i0*s0space;
7948 neighbour_node_number = i0*s1space;
7953 neighbour_node_number = (neigh_n_p-1)*s0space + i0*s1space;
7958 neighbour_node_number = (neigh_n_p-1)*s2space + i0*s1space;
7963 neighbour_node_number = (neigh_n_p-1)*s0space + (neigh_n_p-1)*s2space + i0*s1space;
7969 neighbour_node_number = i0*s2space;
7974 neighbour_node_number = (neigh_n_p-1)*s0space + i0*s2space;
7979 neighbour_node_number = (neigh_n_p-1)*s1space + i0*s2space;
7984 neighbour_node_number = (neigh_n_p-1)*s0space + (neigh_n_p-1)*s1space + i0*s2space;
7991 master_is_not_edge =
true;
7994 if(master_is_not_edge)
break;
7997 master_node_number.push_back(neighbour_node_number);
7998 master_node_pt.push_back(neighbour_node_pt);
8004 if(master_is_not_edge)
8007 for(
unsigned i0=0; i0<neigh_n_p; i0++)
8010 for(
unsigned i1=0; i1<neigh_n_p; i1++)
8012 const unsigned s0space = 1;
8013 const unsigned s1space = neigh_n_p;
8014 const unsigned s2space = neigh_n_p*neigh_n_p;
8020 neighbour_node_number = i0*s0space + i1*s1space;
8025 neighbour_node_number = (neigh_n_p-1)*s2space + i0*s0space + i1*s1space;
8030 neighbour_node_number = i0*s0space + i1*s2space;
8035 neighbour_node_number = (neigh_n_p-1)*s1space + i0*s0space + i1*s2space;
8040 neighbour_node_number = i0*s1space + i1*s2space;
8045 neighbour_node_number = (neigh_n_p-1)*s0space + i0*s1space + i1*s2space;
8052 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8053 OOMPH_EXCEPTION_LOCATION);
8057 master_node_number.push_back(neighbour_node_number);
8058 master_node_pt.push_back(neighbour_node_pt);
8064 unsigned local_node_number=0;
8065 Node* local_node_pt=0;
8068 for(
unsigned i0=0; i0<my_n_p; i0++)
8073 const unsigned s0space = 1;
8074 const unsigned s1space = my_n_p;
8075 const unsigned s2space = my_n_p*my_n_p;
8083 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
8084 s_fraction[1] = 0.0;
8085 s_fraction[2] = 0.0;
8086 local_node_number = i0*s0space;
8088 = this->interpolating_node_pt(local_node_number,value_id);
8092 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
8093 s_fraction[1] = 1.0;
8094 s_fraction[2] = 0.0;
8095 local_node_number = (my_n_p-1)*s1space + i0*s0space;
8097 = this->interpolating_node_pt(local_node_number,value_id);
8101 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
8102 s_fraction[1] = 0.0;
8103 s_fraction[2] = 1.0;
8104 local_node_number = (my_n_p-1)*s2space + i0*s0space;
8106 = this->interpolating_node_pt(local_node_number,value_id);
8110 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
8111 s_fraction[1] = 1.0;
8112 s_fraction[2] = 1.0;
8113 local_node_number = (my_n_p-1)*s1space + (my_n_p-1)*s2space + i0*s0space;
8115 = this->interpolating_node_pt(local_node_number,value_id);
8119 s_fraction[0] = 0.0;
8121 local_one_d_fraction_of_interpolating_node(i0,1,value_id);
8122 s_fraction[2] = 0.0;
8123 local_node_number = i0*s1space;
8125 = this->interpolating_node_pt(local_node_number,value_id);
8128 s_fraction[0] = 1.0;
8130 local_one_d_fraction_of_interpolating_node(i0,1,value_id);
8131 s_fraction[2] = 0.0;
8132 local_node_number = (my_n_p-1)*s0space + i0*s1space;
8134 = this->interpolating_node_pt(local_node_number,value_id);
8137 s_fraction[0] = 0.0;
8139 local_one_d_fraction_of_interpolating_node(i0,1,value_id);
8140 s_fraction[2] = 1.0;
8141 local_node_number = (my_n_p-1)*s2space + i0*s1space;
8143 = this->interpolating_node_pt(local_node_number,value_id);
8146 s_fraction[0] = 1.0;
8148 local_one_d_fraction_of_interpolating_node(i0,1,value_id);
8149 s_fraction[2] = 1.0;
8150 local_node_number = (my_n_p-1)*s0space + (my_n_p-1)*s2space + i0*s1space;
8152 = this->interpolating_node_pt(local_node_number,value_id);
8156 s_fraction[0] = 0.0;
8157 s_fraction[1] = 0.0;
8159 local_one_d_fraction_of_interpolating_node(i0,2,value_id);
8160 local_node_number = i0*s2space;
8162 = this->interpolating_node_pt(local_node_number,value_id);
8165 s_fraction[0] = 1.0;
8166 s_fraction[1] = 0.0;
8168 local_one_d_fraction_of_interpolating_node(i0,2,value_id);
8169 local_node_number = (my_n_p-1)*s0space + i0*s2space;
8171 = this->interpolating_node_pt(local_node_number,value_id);
8174 s_fraction[0] = 0.0;
8175 s_fraction[1] = 1.0;
8177 local_one_d_fraction_of_interpolating_node(i0,2,value_id);
8178 local_node_number = (my_n_p-1)*s1space + i0*s2space;
8180 = this->interpolating_node_pt(local_node_number,value_id);
8183 s_fraction[0] = 1.0;
8184 s_fraction[1] = 1.0;
8186 local_one_d_fraction_of_interpolating_node(i0,2,value_id);
8187 local_node_number = (my_n_p-1)*s0space + (my_n_p-1)*s1space + i0*s2space;
8189 = this->interpolating_node_pt(local_node_number,value_id);
8194 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8195 OOMPH_EXCEPTION_LOCATION);
8199 slave_node_number.push_back(local_node_number);
8200 slave_node_pt.push_back(local_node_pt);
8203 slave_node_s_fraction.push_back(s_fraction);
8207 const unsigned n_slave_nodes = slave_node_pt.size();
8208 const unsigned n_master_nodes = master_node_pt.size();
8209 const unsigned slave_element_nnode_1d = my_n_p;
8210 const unsigned master_element_nnode_1d = neigh_n_p;
8220 slave_nodal_position,slave_weight);
8224 master_nodal_position,master_weight);
8236 const unsigned n_mortar_vertices = 2;
8239 vertex_pos[1] = this->ninterpolating_node_1d(value_id)-1;
8241 for(
unsigned i=0;
i<my_n_p-n_mortar_vertices;
i++)
8242 {non_vertex_pos[
i] =
i+1;}
8248 std::vector<bool> mortar_vertex_on_master_edge(n_mortar_vertices);
8249 for (
unsigned v=0; v<n_mortar_vertices; v++)
8252 mortar_vertex_on_master_edge[v] =
true;
8253 unsigned non_extreme_coordinate = 0;
8255 for(
unsigned i=0;
i<3;
i++)
8258 s_in_neigh[
i] = edge_s_lo_neigh[
i] +
8259 slave_node_s_fraction[vertex_pos[v]][edge_translate_s[
i]]*
8260 (edge_s_hi_neigh[
i] - edge_s_lo_neigh[
i]);
8263 if(std::fabs(std::fabs(s_in_neigh[
i])-1.0) > 1.0e-14)
8265 non_extreme_coordinate++;
8266 if(non_extreme_coordinate>1)
8268 mortar_vertex_on_master_edge[v] =
false;
8276 bool my_edge_coincides_with_master =
true;
8277 for(
unsigned v=0; v<n_mortar_vertices; v++)
8279 my_edge_coincides_with_master = my_edge_coincides_with_master
8280 && mortar_vertex_on_master_edge[v];
8286 for (
unsigned v=0; v<n_mortar_vertices; v++)
8291 if((!my_edge_coincides_with_master) && mortar_vertex_on_master_edge[v])
8295 bool node_is_shared=
false;
8296 for(
unsigned i=0;
i<master_node_pt.size();
i++)
8298 if(slave_node_pt[vertex_pos[v]]==master_node_pt[
i])
8300 node_is_shared=
true;
8314 for(
unsigned i=0;
i<3;
i++)
8316 s_in_neigh[
i] = edge_s_lo_neigh[
i] +
8317 slave_node_s_fraction[vertex_pos[v]][edge_translate_s[
i]]*
8318 (edge_s_hi_neigh[
i] - edge_s_lo_neigh[
i]);
8326 slave_node_pt[vertex_pos[v]]->set_nonhanging();
8330 for(
unsigned m=0; m<n_master_nodes; m++)
8333 if(std::fabs(master_shapes[master_node_number[m]]) < 1.0
e-14)
8336 master_node_zero_weight.push_back(m);
8342 new HangInfo(n_master_nodes-master_node_zero_weight.size());
8343 unsigned mindex = 0;
8344 for(
unsigned m=0; m<n_master_nodes; m++)
8348 for(
unsigned i=0;
i<master_node_zero_weight.size();
i++)
8350 if(m==master_node_zero_weight[
i]) skip =
true;
8355 explicit_hang_pt->
set_master_node_pt(mindex++,master_node_pt[m],master_shapes[master_node_number[m]]);
8366 slave_node_pt[vertex_pos[v]]->set_hanging_pt(explicit_hang_pt,-1);
8385 for(
unsigned m=0; m<explicit_hang_pt->
nmaster(); m++)
8391 "Master has zero weight!",
8392 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8393 OOMPH_EXCEPTION_LOCATION);
8401 if(n_slave_nodes-n_mortar_vertices > 0)
8409 for (
unsigned i=0;
i<shared_node_M.size();
i++)
8410 {shared_node_M[
i].resize(n_slave_nodes-n_mortar_vertices);}
8413 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
8416 psi[
i] = pow(-1.0,
int((slave_element_nnode_1d-1)-
i-1))
8418 slave_nodal_position[non_vertex_pos[
i]]);
8420 diag_M[
i] = psi[
i]*slave_weight[non_vertex_pos[
i]];
8424 for(
unsigned v=0; v<shared_node_M.size(); v++)
8426 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
8429 if (std::fabs(slave_nodal_position[non_vertex_pos[
i]]
8430 - slave_nodal_position[vertex_pos[v]]) >= 1.0e-8 )
8433 psi[
i] = pow(-1.0,
int((slave_element_nnode_1d-1)-
i-1))
8435 slave_nodal_position[vertex_pos[v]])
8436 / (slave_nodal_position[non_vertex_pos[
i]]
8437 - slave_nodal_position[vertex_pos[v]]);
8441 slave_nodal_position[vertex_pos[v]]))
8445 psi[
i] = pow(-1.0,
int((slave_element_nnode_1d-1)-
i-1))
8447 slave_nodal_position[non_vertex_pos[
i]]);
8454 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
8455 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8456 OOMPH_EXCEPTION_LOCATION);
8459 shared_node_M[v][
i] = psi[
i]*slave_weight[vertex_pos[v]];
8468 for(
unsigned i=0;
i<P.size();
i++) {P[
i].resize(n_master_nodes,0.0);}
8482 if (slave_element_nnode_1d >= master_element_nnode_1d)
8485 quadrature_knot = &slave_nodal_position;
8486 quadrature_weight = &slave_weight;
8491 quadrature_knot = &master_nodal_position;
8492 quadrature_weight = &master_weight;
8496 for (
unsigned q=0; q<(*quadrature_weight).size(); q++)
8501 s_on_mortar[0] = (*quadrature_knot)[q];
8504 for(
unsigned k=0; k<n_slave_nodes-n_mortar_vertices; k++)
8507 if (std::fabs(slave_nodal_position[non_vertex_pos[k]]-s_on_mortar[0]) >= 1.0
e-08)
8510 psi[k] = pow(-1.0,
int((slave_element_nnode_1d-1)-k-1))
8512 / (slave_nodal_position[non_vertex_pos[k]]-s_on_mortar[0]);
8519 psi[k] = pow(-1.0,
int((slave_element_nnode_1d-1)-k-1))
8527 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
8528 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8529 OOMPH_EXCEPTION_LOCATION);
8535 for(
unsigned i=0;
i<3;
i++)
8537 s_fraction[
i] = (
i==active_coord_index)
8538 ? 0.5*(s_on_mortar[0]+1.0)
8539 : slave_node_s_fraction[vertex_pos[0]][
i];
8544 for(
unsigned i=0;
i<3;
i++)
8546 s_in_neigh[
i] = edge_s_lo_neigh[
i] + s_fraction[edge_translate_s[
i]]*
8547 (edge_s_hi_neigh[
i] - edge_s_lo_neigh[
i]);
8554 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
8556 for(
unsigned j=0; j<n_master_nodes; j++)
8558 P[
i][j] += master_shapes[master_node_number[j]]*psi[
i]*(*quadrature_weight)[q];
8580 for(
unsigned v=0; v<n_mortar_vertices; v++)
8584 for(
unsigned i=0;
i<3;
i++)
8586 s_fraction[
i] = (
i==active_coord_index)
8587 ? 0.5*(slave_nodal_position[vertex_pos[v]]+1.0)
8588 : slave_node_s_fraction[vertex_pos[0]][
i];
8593 for(
unsigned i=0;
i<3;
i++)
8595 s_in_neigh[
i] = edge_s_lo_neigh[
i] + s_fraction[edge_translate_s[
i]]*
8596 (edge_s_hi_neigh[
i] - edge_s_lo_neigh[
i]);
8602 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
8604 for(
unsigned k=0; k<n_master_nodes; k++)
8606 P[
i][k] -= master_shapes[master_node_number[k]]*shared_node_M[v][
i];
8625 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
8627 for(
unsigned j=0; j<n_master_nodes; j++)
8648 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
8652 for(
unsigned m=0; m<n_master_nodes; m++)
8655 if(std::fabs(P[
i][m]) < 1.0
e-14)
8658 master_node_zero_weight.push_back(m);
8664 new HangInfo(n_master_nodes-master_node_zero_weight.size());
8665 unsigned mindex = 0;
8666 for(
unsigned m=0; m<n_master_nodes; m++)
8670 for(
unsigned k=0; k<master_node_zero_weight.size(); k++)
8672 if(m==master_node_zero_weight[k]) skip =
true;
8677 hang_info_pt[
i]->set_master_node_pt(mindex++,master_node_pt[m],P[
i][m]);
8701 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
8707 bool node_is_really_shared =
false;
8708 for (
unsigned m=0; m<hang_info_pt[
i]->nmaster(); m++)
8711 if (hang_info_pt[
i]->master_node_pt(m)==slave_node_pt[non_vertex_pos[
i]])
8713 node_is_really_shared =
true;
8717 if(std::fabs(hang_info_pt[
i]->master_weight(m)-1.0) > 1.0e-06)
8721 "Something fishy here -- with shared node at a mortar vertex",
8722 "PRefineableQElemen<2,INITIAL_NNODE_1D>t::quad_hang_helper()",
8723 OOMPH_EXCEPTION_LOCATION);
8730 if(!node_is_really_shared)
8732 slave_node_pt[non_vertex_pos[
i]]->set_hanging_pt(hang_info_pt[
i],-1);
8762 if (h_type_slave || p_type_slave)
8766 if(my_face==
B || my_face==
F)
8768 active_coord_index[0] = 0;
8769 active_coord_index[1] = 1;
8771 else if(my_face==
D || my_face==
U)
8773 active_coord_index[0] = 0;
8774 active_coord_index[1] = 2;
8776 else if(my_face==
L || my_face==
R)
8778 active_coord_index[0] = 1;
8779 active_coord_index[1] = 2;
8784 "Cannot transform coordinates",
8785 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8786 OOMPH_EXCEPTION_LOCATION);
8801 const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
8805 unsigned neighbour_node_number=0;
8806 Node* neighbour_node_pt=0;
8809 for(
unsigned i0=0; i0<neigh_n_p; i0++)
8811 for(
unsigned i1=0; i1<neigh_n_p; i1++)
8813 const unsigned s0space = 1;
8814 const unsigned s1space = neigh_n_p;
8815 const unsigned s2space = neigh_n_p*neigh_n_p;
8821 neighbour_node_number = i0*s0space + i1*s1space;
8827 neighbour_node_number = (neigh_n_p-1)*s2space + i0*s0space + i1*s1space;
8833 neighbour_node_number = i0*s0space + i1*s2space;
8839 neighbour_node_number = (neigh_n_p-1)*s1space + i0*s0space + i1*s2space;
8845 neighbour_node_number = i0*s1space + i1*s2space;
8851 neighbour_node_number = (neigh_n_p-1)*s0space + i0*s1space + i1*s2space;
8858 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8859 OOMPH_EXCEPTION_LOCATION);
8863 master_node_number.push_back(neighbour_node_number);
8864 master_node_pt.push_back(neighbour_node_pt);
8869 unsigned local_node_number=0;
8870 Node* local_node_pt=0;
8873 for(
unsigned i0=0; i0<my_n_p; i0++)
8875 for(
unsigned i1=0; i1<my_n_p; i1++)
8880 const unsigned s0space = 1;
8881 const unsigned s1space = my_n_p;
8882 const unsigned s2space = my_n_p*my_n_p;
8890 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
8892 local_one_d_fraction_of_interpolating_node(i1,1,value_id);
8893 s_fraction[2] = 0.0;
8894 local_node_number = i0*s0space + i1*s1space;
8896 = this->interpolating_node_pt(local_node_number,value_id);
8901 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
8903 local_one_d_fraction_of_interpolating_node(i1,1,value_id);
8904 s_fraction[2] = 1.0;
8905 local_node_number = (my_n_p-1)*s2space + i0*s0space + i1*s1space;
8907 = this->interpolating_node_pt(local_node_number,value_id);
8912 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
8913 s_fraction[1] = 0.0;
8915 local_one_d_fraction_of_interpolating_node(i1,2,value_id);
8916 local_node_number = i0*s0space + i1*s2space;
8918 = this->interpolating_node_pt(local_node_number,value_id);
8923 local_one_d_fraction_of_interpolating_node(i0,0,value_id);
8924 s_fraction[1] = 1.0;
8926 local_one_d_fraction_of_interpolating_node(i1,2,value_id);
8927 local_node_number = (my_n_p-1)*s1space + i0*s0space + i1*s2space;
8929 = this->interpolating_node_pt(local_node_number,value_id);
8933 s_fraction[0] = 0.0;
8935 local_one_d_fraction_of_interpolating_node(i0,1,value_id);
8937 local_one_d_fraction_of_interpolating_node(i1,2,value_id);
8938 local_node_number = i0*s1space + i1*s2space;
8940 = this->interpolating_node_pt(local_node_number,value_id);
8944 s_fraction[0] = 1.0;
8946 local_one_d_fraction_of_interpolating_node(i0,1,value_id);
8948 local_one_d_fraction_of_interpolating_node(i1,2,value_id);
8949 local_node_number = (my_n_p-1)*s0space + i0*s1space + i1*s2space;
8951 = this->interpolating_node_pt(local_node_number,value_id);
8956 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8957 OOMPH_EXCEPTION_LOCATION);
8961 slave_node_number.push_back(local_node_number);
8962 slave_node_pt.push_back(local_node_pt);
8965 slave_node_s_fraction.push_back(s_fraction);
8970 const unsigned n_slave_nodes = slave_node_pt.size();
8971 const unsigned n_master_nodes = master_node_pt.size();
8972 const unsigned slave_element_nnode_1d = my_n_p;
8973 const unsigned master_element_nnode_1d = neigh_n_p;
9004 slave_nodal_position_1d,slave_weight_1d);
9008 master_nodal_position_1d,master_weight_1d);
9012 slave_nodal_position(slave_element_nnode_1d*slave_element_nnode_1d);
9013 for(
unsigned i=0;
i<slave_nodal_position.size();
i++)
9014 {slave_nodal_position[
i].resize(2);}
9015 Vector<double> slave_weight(slave_element_nnode_1d*slave_element_nnode_1d);
9017 master_nodal_position(master_element_nnode_1d*master_element_nnode_1d);
9018 for(
unsigned i=0;
i<master_nodal_position.size();
i++)
9019 {master_nodal_position[
i].resize(2);}
9020 Vector<double> master_weight(master_element_nnode_1d*master_element_nnode_1d);
9023 unsigned slave_index = 0;
9024 for (
unsigned i=0;
i<slave_element_nnode_1d;
i++)
9026 for (
unsigned j=0; j<slave_element_nnode_1d; j++)
9028 slave_nodal_position[slave_index][0] = slave_nodal_position_1d[
i];
9029 slave_nodal_position[slave_index][1] = slave_nodal_position_1d[j];
9030 slave_weight[slave_index] = slave_weight_1d[
i]*slave_weight_1d[j];
9034 unsigned master_index = 0;
9035 for (
unsigned i=0;
i<master_element_nnode_1d;
i++)
9037 for (
unsigned j=0; j<master_element_nnode_1d; j++)
9039 master_nodal_position[master_index][0] = master_nodal_position_1d[
i];
9040 master_nodal_position[master_index][1] = master_nodal_position_1d[j];
9041 master_weight[master_index] = master_weight_1d[
i]*master_weight_1d[j];
9061 Vector<unsigned> non_vertex_pos((slave_element_nnode_1d-2)*(slave_element_nnode_1d-2));
9062 unsigned nvindex = 0;
9063 for(
unsigned i=1;
i<slave_element_nnode_1d-1;
i++)
9065 for(
unsigned j=1; j<slave_element_nnode_1d-1; j++)
9067 non_vertex_pos[nvindex++] =
i*slave_element_nnode_1d + j;
9071 for(
unsigned i=0;
i<n_slave_nodes;
i++)
9074 bool node_is_vertex =
true;
9075 for(
unsigned j=0; j<non_vertex_pos.size(); j++)
9077 if(
i==non_vertex_pos[j])
9080 node_is_vertex =
false;
9087 vertex_pos.push_back(
i);
9091 const unsigned n_mortar_vertices = vertex_pos.size();
9095 if(n_slave_nodes-n_mortar_vertices > 0)
9103 for (
unsigned i=0;
i<shared_node_M.size();
i++)
9104 {shared_node_M[
i].resize(n_slave_nodes-n_mortar_vertices);}
9107 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
9114 for(
unsigned dir=0; dir<2; dir++)
9116 unsigned index1d = (dir==0) ?
i :
i % (slave_element_nnode_1d-2);
9118 psi[
i] *= pow(-1.0,
int((slave_element_nnode_1d-1)-index1d-1))
9120 slave_nodal_position[non_vertex_pos[
i]][dir]);
9123 diag_M[
i] = psi[
i]*slave_weight[non_vertex_pos[
i]];
9135 for(
unsigned v=0; v<shared_node_M.size(); v++)
9137 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
9144 for(
unsigned dir=0; dir<2; dir++)
9146 unsigned index1d = (dir==0) ?
i :
i % (slave_element_nnode_1d-2);
9148 if (std::fabs(slave_nodal_position[non_vertex_pos[
i]][dir]
9149 - slave_nodal_position[vertex_pos[v]][dir]) >= 1.0e-8 )
9152 psi[
i] *= pow(-1.0,
int((slave_element_nnode_1d-1)-index1d-1))
9154 slave_nodal_position[vertex_pos[v]][dir])
9155 / (slave_nodal_position[non_vertex_pos[
i]][dir]
9156 - slave_nodal_position[vertex_pos[v]][dir]);
9161 slave_nodal_position[vertex_pos[v]][dir]))
9165 psi[
i] *= pow(-1.0,
int((slave_element_nnode_1d-1)-index1d-1))
9167 slave_nodal_position[non_vertex_pos[
i]][dir]);
9174 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
9175 "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
9176 OOMPH_EXCEPTION_LOCATION);
9180 shared_node_M[v][
i] = psi[
i]*slave_weight[vertex_pos[v]];
9200 for(
unsigned i=0;
i<P.size();
i++) {P[
i].resize(n_master_nodes,0.0);}
9215 if (slave_element_nnode_1d >= master_element_nnode_1d)
9218 quadrature_knot = &slave_nodal_position;
9219 quadrature_weight = &slave_weight;
9224 quadrature_knot = &master_nodal_position;
9225 quadrature_weight = &master_weight;
9229 for (
unsigned q=0; q<(*quadrature_weight).size(); q++)
9233 for(
unsigned i=0;
i<2;
i++)
9234 {s_on_mortar[
i] = (*quadrature_knot)[q][
i];}
9237 for(
unsigned k=0; k<n_slave_nodes-n_mortar_vertices; k++)
9244 for(
unsigned dir=0; dir<2; dir++)
9246 unsigned index1d = (dir==0) ? k : k % (slave_element_nnode_1d-2);
9249 slave_nodal_position[non_vertex_pos[k]][dir]-s_on_mortar[dir])
9253 psi[k] *= pow(-1.0,
int((slave_element_nnode_1d-1)-index1d-1))
9255 slave_element_nnode_1d-1,s_on_mortar[dir])
9256 / (slave_nodal_position[non_vertex_pos[k]][dir]-s_on_mortar[dir]);
9261 slave_element_nnode_1d-1,s_on_mortar[dir]))
9265 psi[k] *= pow(-1.0,
int((slave_element_nnode_1d-1)-index1d-1))
9267 slave_element_nnode_1d-1,s_on_mortar[dir]);
9274 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
9275 "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
9276 OOMPH_EXCEPTION_LOCATION);
9283 for(
unsigned i=0;
i<3;
i++)
9285 if(
i==active_coord_index[0])
9286 {s_fraction[
i] = 0.5*(s_on_mortar[0]+1.0);}
9287 else if(
i==active_coord_index[1])
9288 {s_fraction[
i] = 0.5*(s_on_mortar[1]+1.0);}
9290 {s_fraction[
i] = slave_node_s_fraction[vertex_pos[0]][
i];}
9295 for(
unsigned i=0;
i<3;
i++)
9297 s_in_neigh[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
9298 (s_hi_neigh[
i] - s_lo_neigh[
i]);
9305 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
9307 for(
unsigned j=0; j<n_master_nodes; j++)
9309 P[
i][j] += master_shapes[master_node_number[j]]*psi[
i]*(*quadrature_weight)[q];
9331 for(
unsigned v=0; v<n_mortar_vertices; v++)
9335 for(
unsigned i=0;
i<3;
i++)
9337 if(
i==active_coord_index[0])
9338 {s_fraction[
i] = 0.5*(slave_nodal_position[vertex_pos[v]][0]+1.0);}
9339 else if(
i==active_coord_index[1])
9340 {s_fraction[
i] = 0.5*(slave_nodal_position[vertex_pos[v]][1]+1.0);}
9342 {s_fraction[
i] = slave_node_s_fraction[vertex_pos[0]][
i];}
9347 for(
unsigned i=0;
i<3;
i++)
9349 s_in_neigh[
i] = s_lo_neigh[
i] + s_fraction[translate_s[
i]]*
9350 (s_hi_neigh[
i] - s_lo_neigh[
i]);
9356 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
9358 for(
unsigned k=0; k<n_master_nodes; k++)
9360 P[
i][k] -= master_shapes[master_node_number[k]]*shared_node_M[v][
i];
9378 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
9380 for(
unsigned j=0; j<n_master_nodes; j++)
9400 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
9402 hang_info_pt[
i] =
new HangInfo(n_master_nodes);
9407 for(
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
9409 for(
unsigned j=0; j<n_master_nodes; j++)
9411 hang_info_pt[
i]->set_master_node_pt(j,master_node_pt[j],P[
i][j]);
9417 for (
unsigned i=0;
i<n_slave_nodes-n_mortar_vertices;
i++)
9423 bool node_is_really_shared =
false;
9424 for (
unsigned m=0; m<hang_info_pt[
i]->nmaster(); m++)
9427 if (hang_info_pt[
i]->master_node_pt(m)==slave_node_pt[non_vertex_pos[
i]])
9429 node_is_really_shared =
true;
9433 if(std::fabs(hang_info_pt[
i]->master_weight(m)-1.0) > 1.0e-06)
9437 "Something fishy here -- with shared node at a mortar vertex",
9438 "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
9439 OOMPH_EXCEPTION_LOCATION);
9446 if(!node_is_really_shared)
9448 slave_node_pt[non_vertex_pos[
i]]->set_hanging_pt(hang_info_pt[
i],-1);
9456 for (
unsigned i=0;
i<n_slave_nodes;
i++)
9459 if (slave_node_pt[
i]->is_hanging())
9462 double weight_sum = 0.0;
9463 bool zero_weight =
false;
9464 for (
unsigned m=0; m<slave_node_pt[
i]->hanging_pt()->nmaster(); m++)
9466 weight_sum += slave_node_pt[
i]->hanging_pt()->master_weight(m);
9467 if(std::fabs(slave_node_pt[
i]->hanging_pt()->master_weight(m)) < 1.0
e-14)
9470 oomph_info <<
"In the hanging scheme for slave node " <<
i <<
", master node " << m <<
" has weight " << slave_node_pt[
i]->hanging_pt()->master_weight(m) <<
" < 1.0e-14" << std::endl;
9475 if (std::fabs(weight_sum-1.0) > 1.0
e-08)
9477 oomph_info <<
"Sum of master weights: " << weight_sum << std::endl;
9479 "Weights in hanging scheme do not sum to 1",
9480 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9481 OOMPH_EXCEPTION_LOCATION);
9486 "Zero weights present in hanging schemes",
9487 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9488 OOMPH_EXCEPTION_LOCATION);
9493 for (
unsigned m=0; m<slave_node_pt[
i]->hanging_pt()->nmaster(); m++)
9495 if(slave_node_pt[
i] == slave_node_pt[
i]->hanging_pt()->master_node_pt(m))
9499 "Slave node has itself as a master!",
9500 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9501 OOMPH_EXCEPTION_LOCATION);
9503 if(slave_node_pt[
i]->hanging_pt()->master_node_pt(m)->is_hanging())
9506 Node* new_nod_pt = slave_node_pt[
i]->hanging_pt()->master_node_pt(m);
9511 std::cout <<
"++++++++++++++++++++++++++++++++++++++++" << std::endl;
9512 std::cout <<
" Slave node: " 9513 << slave_node_pt[
i] <<
" = " 9514 << slave_node_pt[
i]->x(0) <<
" " 9515 << slave_node_pt[
i]->x(1) <<
" " 9516 << slave_node_pt[
i]->x(2) <<
" " << std::endl;
9517 std::cout <<
" Master node: " 9518 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m) <<
" = " 9519 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m)->x(0) <<
" " 9520 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m)->x(1) <<
" " 9521 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m)->x(2) <<
" " << std::endl;
9522 std::cout <<
"Master's master node: " 9523 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m)->hanging_pt()->master_node_pt(mm) <<
" = " 9524 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m)->hanging_pt()->master_node_pt(mm)->x(0) <<
" " 9525 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m)->hanging_pt()->master_node_pt(mm)->x(1) <<
" " 9526 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m)->hanging_pt()->master_node_pt(mm)->x(2) <<
" " << std::endl;
9531 std::cout <<
"Hanging node: " 9532 << slave_node_pt[
i]->x(0) <<
" " 9533 << slave_node_pt[
i]->x(1) <<
" " 9534 << slave_node_pt[
i]->x(2) <<
" " 9536 for(
unsigned m_tmp=0; m_tmp<slave_node_pt[
i]->hanging_pt()->nmaster(); m_tmp++)
9538 std::cout <<
" m = " << m_tmp <<
": " 9539 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m_tmp)->x(0) <<
" " 9540 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m_tmp)->x(1) <<
" " 9541 << slave_node_pt[
i]->hanging_pt()->master_node_pt(m_tmp)->x(2) <<
" " 9542 <<
"w = " << slave_node_pt[
i]->hanging_pt()->master_weight(m_tmp) <<
" " 9547 std::cout <<
"Master node " << m <<
" of Hanging node: " 9548 << new_nod_pt->
x(0) <<
" " 9549 << new_nod_pt->
x(1) <<
" " 9550 << new_nod_pt->
x(2) <<
" " 9554 std::cout <<
" mm = " << mm_tmp <<
": " 9564 "Slave node has itself as a master of a master!",
9565 "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9566 OOMPH_EXCEPTION_LOCATION);
9576 bool is_master =
false;
9577 for (
unsigned n=0; n<n_master_nodes; n++)
9579 if (slave_node_pt[
i]==master_node_pt[n])
9611 for(
unsigned i=0;
i<n_slave_nodes;
i++)
9614 if(slave_node_pt[
i]->is_hanging())
9621 this->local_coordinate_of_node(slave_node_number[
i],s_local);
9626 this->interpolated_x(s_local,x_in_neighb);
9630 slave_node_pt[
i]->x(0)=x_in_neighb[0];
9631 slave_node_pt[
i]->x(1)=x_in_neighb[1];
9632 slave_node_pt[
i]->x(2)=x_in_neighb[2];
unsigned nsons() const
Return number of sons (zero if it's a leaf node)
A Generalised Element class.
p-refineable version of RefineableQElement<3,INITIAL_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_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
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 Node * interpolating_node_pt(const unsigned &n, const int &value_id)
In mixed elements, different sets of nodes are used to interpolate different unknowns. This function returns the n-th node that interpolates the value_id-th unknown. Default implementation is that all variables use the positional nodes, i.e. isoparametric elements. Note that any overloaded versions of this function MUST provide a set of nodes for the position, which always has the value_id -1.
double dlegendre(const unsigned &p, const double &x)
Calculates first derivative of Legendre polynomial of degree p at x using three term recursive formul...
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Node * get_node_at_local_coordinate(const Vector< double > &s) const
Return the node at the specified local coordinate.
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.
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
unsigned & p_order()
Access function to P_order.
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...
bool is_neighbour_periodic(const int &direction)
Return whether the neighbour in the particular direction is periodic.
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
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...
unsigned nmaster() const
Return the number of master nodes.
static double nodal_position(const unsigned &n)
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
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.
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 pin(const unsigned &i)
Pin the i-th stored variable.
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...
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Node * get_node_at_local_coordinate(const Vector< double > &s) const
Return the node at the specified local coordinate.
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...
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 pin_position(const unsigned &i)
Pin the nodal position.
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Refineable version of Solid brick elements.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
double ddlegendre(const unsigned &p, const double &x)
Calculates second derivative of Legendre polynomial of degree p at x using three term recursive formu...
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...
Node * get_node_at_local_coordinate(const Vector< double > &s) const
Return the node at the specified local coordinate.
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 nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
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...
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Vector< GeomObject * > & geom_object_pt()
Vector of (pointers to) geometric objects involved in node update function.
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Class that contains data for hanging nodes.
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
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...
static const int OMEGA
Default value for an unassigned neighbour.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Base class for elements that allow MacroElement-based node update.
void local_coordinate_of_node(const unsigned &n, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size.
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...
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 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...
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
p-refineable version of RefineableElement
static void calculate_nodal_positions()
Static function used to populate the stored positions.
p-refineable version of RefineableQElement<2,INITIAL_NNODE_1D>.
Refineable version of Solid quad elements.
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) ...
int son_type() const
Return son type.
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
void get_boundaries(const int &edge, std::set< unsigned > &boundaries) const
unsigned nnode() const
Return the number of nodes.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
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...
TreeRoot *& root_pt()
Return pointer to root of the tree.
Class that returns the shape functions associated with legendre.
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...