hp_refineable_elements.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Non-inline member functions for hp-refineable elements
31 
32 //oomph-lib includes
33 #include "algebraic_elements.h"
35 #include "hp_refineable_elements.h"
36 //#include "shape.h"
37 
38 namespace oomph
39 {
40 
41 ////////////////////////////////////////////////////////////////
42 // 1D PRefineableQElements
43 ////////////////////////////////////////////////////////////////
44 
45 /// Get local coordinates of node j in the element; vector sets its own size
46 template<unsigned INITIAL_NNODE_1D>
48 local_coordinate_of_node(const unsigned& n, Vector<double>& s) const
49  {
50  s.resize(1);
51 
52  switch(this->nnode_1d())
53  {
54  case 2:
57  break;
58  case 3:
61  break;
62  case 4:
65  break;
66  case 5:
69  break;
70  case 6:
73  break;
74  case 7:
77  break;
78  default:
79  oomph_info << "\n ERROR: Exceeded maximum polynomial order for";
80  oomph_info << "\n shape functions." << std::endl;
81  break;
82  }
83  }
84 
85 /// Get the local fractino of node j in the element
86 template<unsigned INITIAL_NNODE_1D>
88 local_fraction_of_node(const unsigned &n, Vector<double> &s_fraction)
89  {
90  this->local_coordinate_of_node(n,s_fraction);
91  s_fraction[0] = 0.5*(s_fraction[0] + 1.0);
92  }
93 
94 template<unsigned INITIAL_NNODE_1D>
96 local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
97  {
98  switch(this->nnode_1d())
99  {
100  case 2:
102  return 0.5*(OneDimensionalLegendreShape<2>::nodal_position(n1d) + 1.0);
103  case 3:
105  return 0.5*(OneDimensionalLegendreShape<3>::nodal_position(n1d) + 1.0);
106  case 4:
108  return 0.5*(OneDimensionalLegendreShape<4>::nodal_position(n1d) + 1.0);
109  case 5:
111  return 0.5*(OneDimensionalLegendreShape<5>::nodal_position(n1d) + 1.0);
112  case 6:
114  return 0.5*(OneDimensionalLegendreShape<6>::nodal_position(n1d) + 1.0);
115  case 7:
117  return 0.5*(OneDimensionalLegendreShape<7>::nodal_position(n1d) + 1.0);
118  default:
119  std::ostringstream error_message;
120  error_message <<"\nERROR: Exceeded maximum polynomial order for";
121  error_message <<"\n shape functions.\n";
122  throw OomphLibError(error_message.str(),
123  OOMPH_CURRENT_FUNCTION,
124  OOMPH_EXCEPTION_LOCATION);
125  return 0.0;
126  }
127  }
128 
129 
130 //==================================================================
131 /// Return the node at the specified local coordinate
132 //==================================================================
133 template<unsigned INITIAL_NNODE_1D>
136 {
137  //Load the tolerance into a local variable
139  //There is one possible index.
140  Vector<int> index(1);
141 
142  // Determine the index
143  // -------------------
144 
145  // If we are at the lower limit, the index is zero
146  if(std::fabs(s[0] + 1.0) < tol)
147  {
148  index[0] = 0;
149  }
150  // If we are at the upper limit, the index is the number of nodes minus 1
151  else if(std::fabs(s[0] - 1.0) < tol)
152  {
153  index[0] = this->nnode_1d()-1;
154  }
155  // Otherwise, we have to calculate the index in general
156  else
157  {
158  // Compute Gauss-Lobatto-Legendre node positions
159  Vector<double> z;
160  Orthpoly::gll_nodes(this->nnode_1d(), z);
161  // Loop over possible internal nodal positions
162  for (unsigned n=1; n<this->nnode_1d()-1; n++)
163  {
164  if (std::fabs(z[n] - s[0]) < tol)
165  {
166  index[0] = n;
167  break;
168  }
169  }
170  // No matching nodes
171  return 0;
172  }
173  // If we've got here we have a node, so let's return a pointer to it
174  return this->node_pt(index[0]);
175 }
176 
177 //===================================================================
178 /// If a neighbouring element's son has already created a node at
179 /// a position corresponding to the local fractional position within the
180 /// present element, s_fraction, return
181 /// a pointer to that node. If not, return NULL (0). If the node is
182 /// periodic the flag is_periodic will be true
183 //===================================================================
184 template<unsigned INITIAL_NNODE_1D>
187  bool &is_periodic)
188 {
189  // Not possible in 1D case, so return null pointer
190  return 0;
191 }
192 
193 //==================================================================
194 /// Set the correct p-order of the element based on that of its
195 /// father. Then construct an integration scheme of the correct order.
196 /// If an adopted father is specified, information from this is
197 /// used instead of using the father found from the tree.
198 //==================================================================
199 template<unsigned INITIAL_NNODE_1D>
201 initial_setup(Tree* const &adopted_father_pt, const unsigned &initial_p_order)
202 {
203  //Storage for pointer to my father (in binarytree impersonation)
204  BinaryTree* father_pt = 0;
205 
206  // Check if an adopted father has been specified
207  if (adopted_father_pt!=0)
208  {
209  //Get pointer to my father (in binarytree impersonation)
210  father_pt = dynamic_cast<BinaryTree*>(adopted_father_pt);
211  }
212  // Check if element is in a tree
213  else if (Tree_pt!=0)
214  {
215  //Get pointer to my father (in binarytree impersonation)
216  father_pt = dynamic_cast<BinaryTree*>(binary_tree_pt()->father_pt());
217  }
218  //else
219  // {
220  // throw OomphLibError(
221  // "Element not in a tree, and no adopted father has been specified!",
222  // OOMPH_CURRENT_FUNCTION,
223  // OOMPH_EXCEPTION_LOCATION);
224  // }
225 
226  // Check if element has father
227  if (father_pt!=0 || initial_p_order!=0)
228  {
229  if(father_pt!=0)
230  {
233  (father_pt->object_pt());
234  if (father_el_pt!=0)
235  {
236  unsigned father_p_order = father_el_pt->p_order();
237  // Set p-order to that of father
238  P_order = father_p_order;
239  }
240  }
241  else
242  {
243  P_order = initial_p_order;
244  }
245 
246  // Now sort out the element...
247  // (has p nodes)
248  unsigned new_n_node = P_order;
249 
250  // Allocate new space for Nodes (at the element level)
251  this->set_n_node(new_n_node);
252 
253  // Set integration scheme
254  delete this->integral_pt();
255  switch(P_order)
256  {
257  case 2:
258  this->set_integration_scheme(new GaussLobattoLegendre<1,2>);
259  break;
260  case 3:
261  this->set_integration_scheme(new GaussLobattoLegendre<1,3>);
262  break;
263  case 4:
264  this->set_integration_scheme(new GaussLobattoLegendre<1,4>);
265  break;
266  case 5:
267  this->set_integration_scheme(new GaussLobattoLegendre<1,5>);
268  break;
269  case 6:
270  this->set_integration_scheme(new GaussLobattoLegendre<1,6>);
271  break;
272  case 7:
273  this->set_integration_scheme(new GaussLobattoLegendre<1,7>);
274  break;
275  default:
276  std::ostringstream error_message;
277  error_message <<"\nERROR: Exceeded maximum polynomial order for";
278  error_message <<"\n integration scheme.\n";
279  throw OomphLibError(error_message.str(),
280  OOMPH_CURRENT_FUNCTION,
281  OOMPH_EXCEPTION_LOCATION);
282  }
283  }
284 }
285 
286 //==================================================================
287 /// Check the father element for any required nodes which
288 /// already exist
289 //==================================================================
290 template<unsigned INITIAL_NNODE_1D>
292  Mesh*& mesh_pt,
293  Vector<Node*>& new_node_pt)
294 {
295  /*
296  //Pointer to my father (in binarytree impersonation)
297  BinaryTree* father_pt = dynamic_cast<BinaryTree*>(binary_tree_pt()->father_pt());
298 
299  // Check if element has father
300  if (father_pt!=0)
301  {
302  PRefineableQElement<1>* father_el_pt =
303  dynamic_cast<PRefineableQElement<1>*>
304  (this->tree_pt()->father_pt()->object_pt());
305  if (father_el_pt!=0)
306  {
307  // Pre-build actions
308  //??
309  }
310  else
311  {
312  std::ostringstream error_message;
313  error_message <<"\nERROR: Dynamic cast failed!\n";
314  throw OomphLibError(error_message.str(),
315  OOMPH_CURRENT_FUNCTION,
316  OOMPH_EXCEPTION_LOCATION);
317  }
318  }
319  */
320 }
321 
322 //==================================================================
323 /// p-refine the element inc times. (If inc<0 then p-unrefinement
324 /// is performed.)
325 //==================================================================
326 template<unsigned INITIAL_NNODE_1D>
328  const int &inc,
329  Mesh* const &mesh_pt,
330  GeneralisedElement* const &clone_pt)
331 {
332  //Cast clone to correct type
334  = dynamic_cast<PRefineableQElement<1,INITIAL_NNODE_1D>*>(clone_pt);
335 
336  //Check if we can use it
337  if(clone_el_pt==0)
338  {
339  throw OomphLibError(
340  "Cloned copy must be a PRefineableQElement<1,INITIAL_NNODE_1D>!",
341  OOMPH_CURRENT_FUNCTION,
342  OOMPH_EXCEPTION_LOCATION);
343  }
344 #ifdef PARANOID
345  //Clone exists, so check that it is infact a copy of me
346  else
347  {
348  //Flag to keep track of check
349  bool clone_is_ok = true;
350 
351  //Does the clone have the correct p-order?
352  clone_is_ok = clone_is_ok && (clone_el_pt->p_order() == this->p_order());
353 
354  if(!clone_is_ok)
355  {
356  std::ostringstream err_stream;
357  err_stream << "Clone element has a different p-order from me,"
358  << std::endl
359  << "but it is supposed to be a copy!" << std::endl;
360 
361  throw OomphLibError(err_stream.str(),
362  OOMPH_CURRENT_FUNCTION,
363  OOMPH_EXCEPTION_LOCATION);
364  }
365 
366  //Does the clone have the same number of nodes as me?
367  clone_is_ok = clone_is_ok && (clone_el_pt->nnode() == this->nnode());
368 
369  if(!clone_is_ok)
370  {
371  std::ostringstream err_stream;
372  err_stream << "Clone element has a different number of nodes from me,"
373  << std::endl
374  << "but it is supposed to be a copy!" << std::endl;
375 
376  throw OomphLibError(err_stream.str(),
377  OOMPH_CURRENT_FUNCTION,
378  OOMPH_EXCEPTION_LOCATION);
379  }
380 
381  //Does the clone have the same nodes as me?
382  for(unsigned n=0; n<this->nnode(); n++)
383  {
384  clone_is_ok = clone_is_ok && (clone_el_pt->node_pt(n) == this->node_pt(n));
385  }
386 
387  if(!clone_is_ok)
388  {
389  std::ostringstream err_stream;
390  err_stream << "The nodes belonging to the clone element are different"
391  << std::endl
392  << "from mine, but it is supposed to be a copy!" << std::endl;
393 
394  throw OomphLibError(err_stream.str(),
395  OOMPH_CURRENT_FUNCTION,
396  OOMPH_EXCEPTION_LOCATION);
397  }
398 
399  //If we get to here then the clone has all the information we require
400  }
401 #endif
402 
403  // Currently we can't handle the case of generalised coordinates
404  // since we haven't established how they should be interpolated.
405  // Buffer this case:
406  if (clone_el_pt->node_pt(0)->nposition_type()!=1)
407  {
408  throw OomphLibError(
409  "Can't handle generalised nodal positions (yet).",
410  OOMPH_CURRENT_FUNCTION,
411  OOMPH_EXCEPTION_LOCATION);
412  }
413 
414  // Timestepper should be the same for all nodes -- use it
415  // to create timesteppers for new nodes
416  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
417 
418  // Get number of history values (incl. present)
419  unsigned ntstorage = time_stepper_pt->ntstorage();
420 
421  // Increment p-order of the element
422  P_order += inc;
423 
424  // Change integration scheme
425  delete this->integral_pt();
426  switch(P_order)
427  {
428  case 2:
429  this->set_integration_scheme(new GaussLobattoLegendre<1,2>);
430  break;
431  case 3:
432  this->set_integration_scheme(new GaussLobattoLegendre<1,3>);
433  break;
434  case 4:
435  this->set_integration_scheme(new GaussLobattoLegendre<1,4>);
436  break;
437  case 5:
438  this->set_integration_scheme(new GaussLobattoLegendre<1,5>);
439  break;
440  case 6:
441  this->set_integration_scheme(new GaussLobattoLegendre<1,6>);
442  break;
443  case 7:
444  this->set_integration_scheme(new GaussLobattoLegendre<1,7>);
445  break;
446  default:
447  std::ostringstream error_message;
448  error_message <<"\nERROR: Exceeded maximum polynomial order for";
449  error_message <<"\n integration scheme.\n";
450  throw OomphLibError(error_message.str(),
451  OOMPH_CURRENT_FUNCTION,
452  OOMPH_EXCEPTION_LOCATION);
453  }
454 
455  // Allocate new space for Nodes (at the element level)
456  this->set_n_node(P_order);
457 
458  // Copy vertex nodes and create new internal nodes
459  //------------------------------------------------
460 
461  // Setup vertex coordinates in element:
462  //-------------------------------------
463  Vector<double> s_left(1);
464  Vector<double> s_right(1);
465  s_left[0]=-1.0;
466  s_right[0]= 1.0;
467 
468  //Local coordinate in element
469  Vector<double> s(1);
470 
471  Vector<double> s_fraction(1);
472  // Loop over all nodes in the element
473  for(unsigned n=0;n<P_order;n++)
474  {
475  // Get the fractional position (in the current element) of the node
476  // in the direction of s[0]
477  s_fraction[0] = this->local_one_d_fraction_of_node(n,0);
478 
479  // Evaluate the local coordinate of the node in the father element
480  s[0] = s_left[0] + (s_right[0]-s_left[0])*s_fraction[0];
481 
482  // Initialise flag: So far, this node hasn't been built or copied yet
483  bool node_done = false;
484 
485  //Get the pointer to the node in this element (or rather, its clone),
486  //returns NULL if there is not node
487  Node* created_node_pt = clone_el_pt->get_node_at_local_coordinate(s);
488 
489  // Does this node already exist in this element?
490  //----------------------------------------------
491  if(created_node_pt!=0)
492  {
493  // Copy node across
494  this->node_pt(n) = created_node_pt;
495 
496  // Make sure that we update the values at the node so that they
497  // are consistent with the present representation. This is only
498  // needed for mixed interpolation where the value at the father
499  // could now become active.
500 
501  // Loop over all history values
502  for(unsigned t=0;t<ntstorage;t++)
503  {
504  // Get values from father element
505  // Note: get_interpolated_values() sets Vector size itself
506  Vector<double> prev_values;
507  clone_el_pt->get_interpolated_values(t,s,prev_values);
508  //Find the minimum number of values
509  //(either those stored at the node, or those returned by
510  // the function)
511  unsigned n_val_at_node = created_node_pt->nvalue();
512  unsigned n_val_from_function = prev_values.size();
513  // Use the ternary conditional operator here
514  unsigned n_var = n_val_at_node < n_val_from_function ?
515  n_val_at_node : n_val_from_function;
516  // Assign the values that we can
517  for(unsigned k=0;k<n_var;k++)
518  {
519  created_node_pt->set_value(t,k,prev_values[k]);
520  }
521  }
522 
523  // Indicate that node has been created by copy
524  node_done = true;
525  }
526 
527  // Node does not exist in this element
528  //------------------------------------
529 
530  // If the node has not been built anywhere ---> build it here
531  if(!node_done)
532  {
533  // In a 1D mesh any node which lies on the boundary must exist in
534  // the initial (coarse) mesh, so any newly-built nodes cannot be
535  // boundary nodes. Therefore we always create a normal "bulk" node.
536 
537  // Create node and set the pointer to it from the element
538  created_node_pt = this->construct_node(n,time_stepper_pt);
539 
540  //Now we set the position and values at the newly created node
541 
542  // Now we set the position and values at the newly created node.
543  // In the first instance use macro element or FE representation
544  // to create past and present nodal positions.
545  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC ELEMENTS AS NOT
546  // ALL OF THEM NECESSARILY IMPLEMENT NONTRIVIAL NODE UPDATE
547  // FUNCTIONS. CALLING THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL
548  // LEAVE THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
549  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL NOT ASSIGN SENSIBLE
550  // INITIAL POSITIONS!)
551 
552  // Loop over history values
553  for(unsigned t=0;t<ntstorage;t++)
554  {
555  // Allocate storage for the previous position of the node
556  Vector<double> x_prev(1);
557 
558  // Get position from father element -- this uses the macro
559  // element representation if appropriate.
560  clone_el_pt->get_x(t,s,x_prev);
561 
562  // Set the previous position of the new node
563  created_node_pt->x(t,0) = x_prev[0];
564 
565  // Allocate storage for the previous values at the node
566  // NOTE: the size of this vector is equal to the number of values
567  // (e.g. 3 velocity components and 1 pressure, say)
568  // associated with each node and NOT the number of history values
569  // which the node stores!
570  Vector<double> prev_values;
571 
572  // Get values from father element
573  // Note: get_interpolated_values() sets Vector size itself.
574  clone_el_pt->get_interpolated_values(t,s,prev_values);
575 
576  // Determine the number of values at the new node
577  const unsigned n_value = created_node_pt->nvalue();
578 
579  // Loop over all values and set the previous values
580  for(unsigned v=0;v<n_value;v++)
581  {
582  created_node_pt->set_value(t,v,prev_values[v]);
583  }
584  } // End of loop over history values
585 
586  // Add new node to mesh (if requested)
587  if(mesh_pt!=0)
588  {
589  mesh_pt->add_node_pt(created_node_pt);
590  }
591 
592  AlgebraicElementBase* alg_el_pt=
593  dynamic_cast<AlgebraicElementBase*>(this);
594 
595  //If we do have an algebraic element
596  if(alg_el_pt!=0)
597  {
598  std::string error_message =
599  "Have not implemented p-refinement for";
600  error_message +=
601  "Algebraic p-refineable elements yet\n";
602 
603  throw
604  OomphLibError(error_message,
605  "PRefineableQElement<1,INITIAL_NNODE_1D>::p_refine()",
606  OOMPH_EXCEPTION_LOCATION);
607  }
608 
609  } // End of case where we build the node ourselves
610 
611  // Check if the element is an algebraic element
612  AlgebraicElementBase* alg_el_pt=
613  dynamic_cast<AlgebraicElementBase*>(this);
614 
615  // If the element is an algebraic element, setup
616  // node position (past and present) from algebraic node update
617  // function. This over-writes previous assingments that
618  // were made based on the macro-element/FE representation.
619  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
620  // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
621  // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
622  // INFO FOR *ALL* ROOT ELEMENTS!
623  if (alg_el_pt!=0)
624  {
625  // Build algebraic node update info for new node
626  // This sets up the node update data for all node update
627  // functions that are shared by all nodes in the father
628  // element
629  alg_el_pt->setup_algebraic_node_update(this->node_pt(n),s,
630  clone_el_pt);
631  }
632 
633  } // End of loop over all nodes in element
634 
635 
636  // If the element is a MacroElementNodeUpdateElement, set
637  // the update parameters for the current element's nodes --
638  // all this needs is the vector of (pointers to the)
639  // geometric objects that affect the MacroElement-based
640  // node update -- this needs to be done to set the node
641  // update info for newly created nodes
642  MacroElementNodeUpdateElementBase* clone_m_el_pt=dynamic_cast<
643  MacroElementNodeUpdateElementBase*>(clone_el_pt);
644  if (clone_m_el_pt!=0)
645  {
646  // Get vector of geometric objects from father (construct vector
647  // via copy operation)
648  Vector<GeomObject*> geom_object_pt(clone_m_el_pt->geom_object_pt());
649 
650  // Cast current element to MacroElementNodeUpdateElement:
651  MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
653 
654 #ifdef PARANOID
655  if (m_el_pt==0)
656  {
657  std::string error_message =
658  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
659  error_message +=
660  "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
661  error_message += "then I should be too....\n";
662 
663  throw OomphLibError(error_message,
664  OOMPH_CURRENT_FUNCTION,
665  OOMPH_EXCEPTION_LOCATION);
666  }
667 #endif
668  // Build update info by passing vector of geometric objects:
669  // This sets the current element to be the update element
670  // for all of the element's nodes -- this is reversed
671  // if the element is ever un-refined in the father element's
672  // rebuild_from_sons() function which overwrites this
673  // assignment to avoid nasty segmentation faults that occur
674  // when a node tries to update itself via an element that no
675  // longer exists...
676  m_el_pt->set_node_update_info(geom_object_pt);
677  }
678 
679 
680  // Loop over all nodes in element again, to re-set the positions
681  // This must be done using the new element's macro-element
682  // representation, rather than the old version which may be
683  // of a different p-order!
684  for(unsigned n=0;n<P_order;n++)
685  {
686  //Get the fractional position of the node in the direction of s[0]
687  s_fraction[0] = this->local_one_d_fraction_of_node(n,0);
688  // Local coordinate
689  s[0] = s_left[0] + (s_right[0]-s_left[0])*s_fraction[0];
690 
691  // Loop over # of history values
692  for (unsigned t=0;t<ntstorage;t++)
693  {
694  // Get position from father element -- this uses the macro
695  // element representation if appropriate. If the node
696  // turns out to be a hanging node later on, then
697  // its position gets adjusted in line with its
698  // hanging node interpolation.
699  Vector<double> x_prev(1);
700  this->get_x(t,s,x_prev);
701 
702  // Set previous positions of the new node
703  this->node_pt(n)->x(t,0) = x_prev[0];
704  }
705  }
706 
707  // Not necessary to delete the old nodes since all original nodes are in the
708  // current mesh and so will be pruned as part of the mesh adaption process.
709 
710  // Do any further-build required
711  this->further_build();
712 }
713 
714 //=======================================================================
715 ///Shape functions for PRefineableQElement<DIM>
716 //=======================================================================
717 template<unsigned INITIAL_NNODE_1D>
719 shape(const Vector<double> &s, Shape &psi) const
720 {
721  switch(p_order())
722  {
723  case 2:
724  {
725  // Calculate nodal positions
727  // Create one-dim shape functions
729  //Now let's loop over the nodal points in the element
730  //and copy the values back in
731  for(unsigned i=0;i<p_order();i++) {psi(i) = psi1[i];}
732  break;
733  }
734  case 3:
735  {
736  // Calculate nodal positions
738  // Create one-dim shape functions
740  //Now let's loop over the nodal points in the element
741  //and copy the values back in
742  for(unsigned i=0;i<p_order();i++) {psi(i) = psi1[i];}
743  break;
744  }
745  case 4:
746  {
747  // Calculate nodal positions
749  // Create one-dim shape functions
751  //Now let's loop over the nodal points in the element
752  //and copy the values back in
753  for(unsigned i=0;i<p_order();i++) {psi(i) = psi1[i];}
754  break;
755  }
756  case 5:
757  {
758  // Calculate nodal positions
760  // Create one-dim shape functions
762  //Now let's loop over the nodal points in the element
763  //and copy the values back in
764  for(unsigned i=0;i<p_order();i++) {psi(i) = psi1[i];}
765  break;
766  }
767  case 6:
768  {
769  // Calculate nodal positions
771  // Create one-dim shape functions
773  //Now let's loop over the nodal points in the element
774  //and copy the values back in
775  for(unsigned i=0;i<p_order();i++) {psi(i) = psi1[i];}
776  break;
777  }
778  case 7:
779  {
780  // Calculate nodal positions
782  // Create one-dim shape functions
784  //Now let's loop over the nodal points in the element
785  //and copy the values back in
786  for(unsigned i=0;i<p_order();i++) {psi(i) = psi1[i];}
787  break;
788  }
789  default:
790  oomph_info << "\n ERROR: PRefineableQElement::shape() exceeded maximum";
791  oomph_info << "\n polynomial order for shape functions." << std::endl;
792  }
793 
794 }
795 
796 //=======================================================================
797 ///Derivatives of shape functions for PRefineableQElement<DIM>
798 //=======================================================================
799 template<unsigned INITIAL_NNODE_1D>
801 dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsi) const
802 {
803  switch(p_order())
804  {
805  case 2:
806  {
807  // Calculate nodal positions
809  //Call the shape functions and derivatives
811  OneDimensionalLegendreDShape<2> dpsi1ds(s[0]);
812  // Loop over shapes and copy across
813  for(unsigned i=0;i<p_order();i++)
814  {
815  psi(i) = psi1[i];
816  dpsi(i,0) = dpsi1ds[i];
817  }
818 
819  break;
820  }
821  case 3:
822  {
823  // Calculate nodal positions
825  //Call the shape functions and derivatives
827  OneDimensionalLegendreDShape<3> dpsi1ds(s[0]);
828  // Loop over shapes and copy across
829  for(unsigned i=0;i<p_order();i++)
830  {
831  psi(i) = psi1[i];
832  dpsi(i,0) = dpsi1ds[i];
833  }
834  break;
835  }
836  case 4:
837  {
838  // Calculate nodal positions
840  //Call the shape functions and derivatives
842  OneDimensionalLegendreDShape<4> dpsi1ds(s[0]);
843  // Loop over shapes and copy across
844  for(unsigned i=0;i<p_order();i++)
845  {
846  psi(i) = psi1[i];
847  dpsi(i,0) = dpsi1ds[i];
848  }
849  break;
850  }
851  case 5:
852  {
853  // Calculate nodal positions
855  //Call the shape functions and derivatives
857  OneDimensionalLegendreDShape<5> dpsi1ds(s[0]);
858  // Loop over shapes and copy across
859  for(unsigned i=0;i<p_order();i++)
860  {
861  psi(i) = psi1[i];
862  dpsi(i,0) = dpsi1ds[i];
863  }
864  break;
865  }
866  case 6:
867  {
868  // Calculate nodal positions
870  //Call the shape functions and derivatives
872  OneDimensionalLegendreDShape<6> dpsi1ds(s[0]);
873  // Loop over shapes and copy across
874  for(unsigned i=0;i<p_order();i++)
875  {
876  psi(i) = psi1[i];
877  dpsi(i,0) = dpsi1ds[i];
878  }
879  break;
880  }
881  case 7:
882  {
883  // Calculate nodal positions
885  //Call the shape functions and derivatives
887  OneDimensionalLegendreDShape<7> dpsi1ds(s[0]);
888  // Loop over shapes and copy across
889  for(unsigned i=0;i<p_order();i++)
890  {
891  psi(i) = psi1[i];
892  dpsi(i,0) = dpsi1ds[i];
893  }
894  break;
895  }
896  default:
897  oomph_info << "\n ERROR: PRefineableQElement::dshape_local() exceeded maximum";
898  oomph_info << "\n polynomial order for shape functions." << std::endl;
899  }
900 
901 }
902 
903 //=======================================================================
904 /// Second derivatives of shape functions for PRefineableQElement<DIM>
905 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
906 //=======================================================================
907 template<unsigned INITIAL_NNODE_1D>
909 d2shape_local(const Vector<double> &s, Shape &psi, DShape &dpsids,
910  DShape &d2psids) const
911 {
912  std::ostringstream error_message;
913  error_message <<"\nd2shape_local currently not implemented for this element\n";
914  throw OomphLibError(error_message.str(),
915  OOMPH_CURRENT_FUNCTION,
916  OOMPH_EXCEPTION_LOCATION);
917 }
918 
919 //=================================================================
920 /// Internal function to set up the hanging nodes on a particular
921 /// edge of the element. (Not required in 1D.)
922 //=================================================================
923 template<unsigned INITIAL_NNODE_1D>
925 binary_hang_helper(const int &value_id, const int &my_edge,
926  std::ofstream& output_hangfile)
927 {}
928 
929 //=======================================================================
930 /// Rebuild the element from nodes found in its sons
931 /// Adjusts its p-order to be the maximum of its sons' p-orders
932 //=======================================================================
933 template<unsigned INITIAL_NNODE_1D>
936 {
937  // Get p-orders of sons
938  unsigned n_sons = this->tree_pt()->nsons();
939  Vector<unsigned> son_p_order(n_sons);
940  unsigned max_son_p_order = 0;
941  for (unsigned ison=0;ison<n_sons;ison++)
942  {
943  PRefineableElement* el_pt = dynamic_cast<PRefineableElement*>(this->tree_pt()->son_pt(ison)->object_pt());
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];
946  }
947 
948  unsigned old_Nnode = this->nnode();
949  unsigned old_P_order = this->p_order();
950  // Set p-order of the element
951  this->p_order() = max_son_p_order;
952 
953  // Change integration scheme
954  delete this->integral_pt();
955  switch(this->p_order())
956  {
957  case 2:
958  this->set_integration_scheme(new GaussLobattoLegendre<1,2>);
959  break;
960  case 3:
961  this->set_integration_scheme(new GaussLobattoLegendre<1,3>);
962  break;
963  case 4:
964  this->set_integration_scheme(new GaussLobattoLegendre<1,4>);
965  break;
966  case 5:
967  this->set_integration_scheme(new GaussLobattoLegendre<1,5>);
968  break;
969  case 6:
970  this->set_integration_scheme(new GaussLobattoLegendre<1,6>);
971  break;
972  case 7:
973  this->set_integration_scheme(new GaussLobattoLegendre<1,7>);
974  break;
975  default:
976  std::ostringstream error_message;
977  error_message <<"\nERROR: Exceeded maximum polynomial order for";
978  error_message <<"\n integration scheme.\n";
979  throw OomphLibError(error_message.str(),
980  OOMPH_CURRENT_FUNCTION,
981  OOMPH_EXCEPTION_LOCATION);
982  }
983 
984  // Back up pointers to old nodes before they are lost
985  Vector<Node*> old_node_pt(old_Nnode);
986  for (unsigned n=0; n<old_Nnode; n++)
987  {
988  old_node_pt[n] = this->node_pt(n);
989  }
990 
991  // Allocate new space for Nodes (at the element level)
992  this->set_n_node(this->p_order());
993 
994  // Copy vertex nodes and create new edge and internal nodes
995  //---------------------------------------------------------
996 
997  // Copy vertex nodes
998  this->node_pt(0) = old_node_pt[0];
999  this->node_pt(this->p_order()-1) = old_node_pt[old_P_order-1];
1000 
1001 
1002  //=============================================================
1003  // Below this line is copied from RefineableQSpectralElement<2>
1004 
1005  // The timestepper should be the same for all nodes and node 0 should
1006  // never be deleted.
1007  if(this->node_pt(0)==0)
1008  {
1009  throw OomphLibError("The vertex node (0) does not exist",
1010  OOMPH_CURRENT_FUNCTION,
1011  OOMPH_EXCEPTION_LOCATION);
1012  }
1013 
1014  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
1015 
1016  // Determine number of history values stored
1017  const unsigned ntstorage = time_stepper_pt->ntstorage();
1018 
1019  // Allocate storage for local coordinates
1020  Vector<double> s_fraction(1), s(1);
1021 
1022  // Determine the number of nodes in the element
1023  const unsigned n_node = this->nnode_1d();
1024 
1025  // Loop over the nodes in the element
1026  for(unsigned n=0;n<n_node;n++)
1027  {
1028  // Get the fractional position of the node in the direction of s[0]
1029  s_fraction[0] = this->local_one_d_fraction_of_node(n,0);
1030 
1031  // Determine the local coordinate in the father element
1032  s[0] = -1.0 + 2.0*s_fraction[0];
1033 
1034  // If the node has not been built
1035  if(this->node_pt(n)==0)
1036  {
1037  // Has the node been created by one of its neighbours?
1038  bool is_periodic = false;
1039  Node* created_node_pt =
1040  this->node_created_by_neighbour(s_fraction,is_periodic);
1041 
1042  // If it has, set the pointer
1043  if(created_node_pt!=0)
1044  {
1045  // If the node is periodic
1046  if(is_periodic)
1047  {
1048  throw OomphLibError(
1049  "Cannot handle periodic nodes yet",
1050  OOMPH_CURRENT_FUNCTION,
1051  OOMPH_EXCEPTION_LOCATION);
1052  }
1053  // Non-periodic case, just set the pointer
1054  else
1055  {
1056  this->node_pt(n) = created_node_pt;
1057  }
1058  }
1059  // Otherwise, we need to build it
1060  else
1061  {
1062  // First, find the son element in which the node should live
1063 
1064  // Find coordinate in the son
1065  Vector<double> s_in_son(1);
1066  using namespace BinaryTreeNames;
1067  int son=-10;
1068  // If s_fraction is between 0 and 0.5, we are in the left son
1069  if(s_fraction[0] < 0.5)
1070  {
1071  son = L;
1072  s_in_son[0] = -1.0 + 4.0*s_fraction[0];
1073  }
1074  // Otherwise we are in the right son
1075  else
1076  {
1077  son = R;
1078  s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
1079  }
1080 
1081  // Get the pointer to the son element in which the new node
1082  // would live
1085  this->tree_pt()->son_pt(son)->object_pt());
1086 
1087  // In 1D we should never be rebuilding an element's vertex nodes
1088  // (since they will never be deleted), so throw an error if we
1089  // appear to be doing so
1090 #ifdef PARANOID
1091  if(n==0 || n==n_node-1)
1092  {
1093  std::string error_message =
1094  "I am trying to rebuild one of the (two) vertex nodes in\n";
1095  error_message +=
1096  "this 1D element. It should not have been possible to delete\n";
1097  error_message +=
1098  "either of these!\n";
1099 
1100  throw OomphLibError(
1101  error_message,
1102  OOMPH_CURRENT_FUNCTION,
1103  OOMPH_EXCEPTION_LOCATION);
1104  }
1105 #endif
1106 
1107  // With this in mind we will always be creating normal "bulk" nodes
1108  this->node_pt(n) = this->construct_node(n,time_stepper_pt);
1109 
1110  // Now we set the position and values at the newly created node
1111 
1112  // In the first instance use macro element or FE representation
1113  // to create past and present nodal positions.
1114  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC ELEMENTS AS NOT
1115  // ALL OF THEM NECESSARILY IMPLEMENT NONTRIVIAL NODE UPDATE
1116  // FUNCTIONS. CALLING THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL
1117  // LEAVE THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
1118  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL NOT ASSIGN SENSIBLE
1119  // INITIAL POSITIONS!)
1120 
1121  // Loop over history values
1122  for(unsigned t=0;t<ntstorage;t++)
1123  {
1124  // Allocate storage for the previous position of the node
1125  Vector<double> x_prev(1);
1126 
1127  // Get position from son element -- this uses the macro element
1128  // representation if appropriate
1129  son_el_pt->get_x(t,s_in_son,x_prev);
1130 
1131  // Set the previous position of the new node
1132  this->node_pt(n)->x(t,0) = x_prev[0];
1133 
1134  // Allocate storage for the previous values at the node
1135  // NOTE: the size of this vector is equal to the number of values
1136  // (e.g. 3 velocity components and 1 pressure, say)
1137  // associated with each node and NOT the number of history values
1138  // which the node stores!
1139  Vector<double> prev_values;
1140 
1141  // Get values from son element
1142  // Note: get_interpolated_values() sets Vector size itself.
1143  son_el_pt->get_interpolated_values(t,s_in_son,prev_values);
1144 
1145  // Determine the number of values at the new node
1146  const unsigned n_value = this->node_pt(n)->nvalue();
1147 
1148  // Loop over all values and set the previous values
1149  for(unsigned v=0;v<n_value;v++)
1150  {
1151  this->node_pt(n)->set_value(t,v,prev_values[v]);
1152  }
1153  } // End of loop over history values
1154 
1155  // Add new node to mesh
1156  mesh_pt->add_node_pt(this->node_pt(n));
1157 
1158  } // End of case where we build the node ourselves
1159 
1160  } // End of if this node has not been built
1161  } // End of loop over nodes in element
1162 
1163  // Check if the element is an algebraic element
1164  // This is done on all nodes in the element after reconstruction
1165  // rather than as the nodes are built
1166  AlgebraicElementBase* alg_el_pt =
1167  dynamic_cast<AlgebraicElementBase*>(this);
1168 
1169  // If so, throw error
1170  if(alg_el_pt!=0)
1171  {
1172  std::string error_message =
1173  "Have not implemented rebuilding from sons for";
1174  error_message +=
1175  "Algebraic p-refineable elements yet\n";
1176 
1177  throw
1178  OomphLibError(error_message,
1179  "PRefineableQElement<1,INITIAL_NNODE_1D>::rebuild_from_sons()",
1180  OOMPH_EXCEPTION_LOCATION);
1181  }
1182 
1183 }
1184 
1185 //=================================================================
1186 /// Check inter-element continuity of
1187 /// - nodal positions
1188 /// - (nodally) interpolated function values
1189 //====================================================================
1190 template<unsigned INITIAL_NNODE_1D>
1192 check_integrity(double& max_error)
1193 {
1195 }
1196 
1197 ////////////////////////////////////////////////////////////////
1198 // 2D PRefineableQElements
1199 ////////////////////////////////////////////////////////////////
1200 
1201 /// Get local coordinates of node j in the element; vector sets its own size
1202 template<unsigned INITIAL_NNODE_1D>
1204 local_coordinate_of_node(const unsigned& n, Vector<double>& s) const
1205  {
1206  s.resize(2);
1207  unsigned Nnode_1d = this->nnode_1d();
1208  unsigned n0 = n%Nnode_1d;
1209  unsigned n1 = n/Nnode_1d;
1210 
1211  switch(Nnode_1d)
1212  {
1213  case 2:
1217  break;
1218  case 3:
1222  break;
1223  case 4:
1227  break;
1228  case 5:
1232  break;
1233  case 6:
1237  break;
1238  case 7:
1242  break;
1243  default:
1244  std::ostringstream error_message;
1245  error_message <<"\nERROR: Exceeded maximum polynomial order for";
1246  error_message <<"\n shape functions.\n";
1247  throw OomphLibError(error_message.str(),
1248  OOMPH_CURRENT_FUNCTION,
1249  OOMPH_EXCEPTION_LOCATION);
1250  break;
1251  }
1252  }
1253 
1254 /// Get the local fractino of node j in the element
1255 template<unsigned INITIAL_NNODE_1D>
1257 local_fraction_of_node(const unsigned &n, Vector<double> &s_fraction)
1258  {
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);
1262  }
1263 
1264 template<unsigned INITIAL_NNODE_1D>
1266 local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
1267  {
1268  switch(this->nnode_1d())
1269  {
1270  case 2:
1272  return 0.5*(OneDimensionalLegendreShape<2>::nodal_position(n1d) + 1.0);
1273  case 3:
1275  return 0.5*(OneDimensionalLegendreShape<3>::nodal_position(n1d) + 1.0);
1276  case 4:
1278  return 0.5*(OneDimensionalLegendreShape<4>::nodal_position(n1d) + 1.0);
1279  case 5:
1281  return 0.5*(OneDimensionalLegendreShape<5>::nodal_position(n1d) + 1.0);
1282  case 6:
1284  return 0.5*(OneDimensionalLegendreShape<6>::nodal_position(n1d) + 1.0);
1285  case 7:
1287  return 0.5*(OneDimensionalLegendreShape<7>::nodal_position(n1d) + 1.0);
1288  default:
1289  std::ostringstream error_message;
1290  error_message <<"\nERROR: Exceeded maximum polynomial order for";
1291  error_message <<"\n shape functions.\n";
1292  throw OomphLibError(error_message.str(),
1293  OOMPH_CURRENT_FUNCTION,
1294  OOMPH_EXCEPTION_LOCATION);
1295  return 0.0;
1296  }
1297  }
1298 
1299 
1300 //==================================================================
1301 /// Return the node at the specified local coordinate
1302 //==================================================================
1303 template<unsigned INITIAL_NNODE_1D>
1306 {
1307  unsigned Nnode_1d = this->nnode_1d();
1308  //Load the tolerance into a local variable
1310  //There are two possible indices.
1311  Vector<int> index(2);
1312 
1313  // Loop over indices
1314  for (unsigned i=0; i<2; i++)
1315  {
1316  // Determine the index
1317  // -------------------
1318 
1319  bool is_found=false;
1320 
1321  // If we are at the lower limit, the index is zero
1322  if(std::fabs(s[i] + 1.0) < tol)
1323  {
1324  index[i] = 0;
1325  is_found=true;
1326  }
1327  // If we are at the upper limit, the index is the number of nodes minus 1
1328  else if(std::fabs(s[i] - 1.0) < tol)
1329  {
1330  index[i] = Nnode_1d-1;
1331  is_found=true;
1332  }
1333  // Otherwise, we have to calculate the index in general
1334  else
1335  {
1336  // Compute Gauss-Lobatto-Legendre node positions
1337  Vector<double> z;
1338  Orthpoly::gll_nodes(Nnode_1d, z);
1339  // Loop over possible internal nodal positions
1340  for (unsigned n=1; n<Nnode_1d-1; n++)
1341  {
1342  if (std::fabs(z[n] - s[i]) < tol)
1343  {
1344  index[i] = n;
1345  is_found=true;
1346  break;
1347  }
1348  }
1349  }
1350 
1351  if (!is_found)
1352  {
1353  // No matching nodes
1354  return 0;
1355  }
1356  }
1357  // If we've got here we have a node, so let's return a pointer to it
1358  return this->node_pt(index[0] + Nnode_1d*index[1]);
1359 }
1360 
1361 //===================================================================
1362 /// If a neighbouring element has already created a node at
1363 /// a position corresponding to the local fractional position within the
1364 /// present element, s_fraction, return
1365 /// a pointer to that node. If not, return NULL (0). If the node is
1366 /// periodic the flag is_periodic will be true
1367 //===================================================================
1368 template<unsigned INITIAL_NNODE_1D>
1370 node_created_by_neighbour(const Vector<double> &s_fraction, bool &is_periodic)
1371 {
1372  using namespace QuadTreeNames;
1373 
1374  //Calculate the edges on which the node lies
1375  Vector<int> edges;
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);}
1380 
1381  //Find the number of edges
1382  unsigned n_size = edges.size();
1383  //If there are no edges, then there is no neighbour, return 0
1384  if(n_size==0) {return 0;}
1385 
1386  Vector<unsigned> translate_s(2);
1387  Vector<double> s_lo_neigh(2);
1388  Vector<double> s_hi_neigh(2);
1389  Vector<double> s(2);
1390 
1391  int neigh_edge, diff_level;
1392  bool in_neighbouring_tree;
1393 
1394  //Loop over the edges
1395  for(unsigned j=0;j<n_size;j++)
1396  {
1397  // Find pointer to neighbouring element along edge
1398  QuadTree* neigh_pt;
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);
1402 
1403  // Neighbour exists
1404  if(neigh_pt!=0)
1405  {
1406  // Have any of its vertex nodes been created yet?
1407  // (Must look in incomplete neighbours because after the
1408  // pre-build they may contain pointers to the required nodes. e.g.
1409  // h-refinement of neighbouring linear and quadratic elements)
1410  bool a_vertex_node_is_built = false;
1411  QElement<2,INITIAL_NNODE_1D>* neigh_obj_pt =
1412  dynamic_cast<QElement<2,INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
1413  if(neigh_obj_pt==0)
1414  {
1415  throw
1416  OomphLibError("Not a quad element!",
1417  "PRefineableQElement<2,INITIAL_NNODE_1D>::node_created_by_neighbour()",
1418  OOMPH_EXCEPTION_LOCATION);
1419  }
1420  for(unsigned vnode=0; vnode<neigh_obj_pt->nvertex_node(); vnode++)
1421  {
1422  if(neigh_obj_pt->vertex_node_pt(vnode)!=0)
1423  a_vertex_node_is_built = true;
1424  break;
1425  }
1426  if(a_vertex_node_is_built)
1427  {
1428  //We now need to translate the nodal location
1429  //as defined in terms of the fractional coordinates of the present
1430  //element into those of its neighbour
1431 
1432  //Calculate the local coordinate in the neighbour
1433  //Note that we need to use the translation scheme in case
1434  //the local coordinates are swapped in the neighbour.
1435  for(unsigned i=0;i<2;i++)
1436  {
1437  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
1438  (s_hi_neigh[i] - s_lo_neigh[i]);
1439  }
1440 
1441  //Find the node in the neighbour
1442  Node* neighbour_node_pt =
1443  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
1444 
1445  //If there is a node, return it
1446  if(neighbour_node_pt!=0)
1447  {
1448  //Now work out whether it's a periodic boundary
1449  //only possible if we have moved into a neighbouring tree
1450  if(in_neighbouring_tree)
1451  {
1452  //Return whether the neighbour is periodic
1453  is_periodic =
1454  quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1455  }
1456  //Return the pointer to the neighbouring node
1457  return neighbour_node_pt;
1458  }
1459  }
1460  }
1461  }
1462  //Node not found, return null
1463  return 0;
1464 }
1465 
1466 //===================================================================
1467 /// If a neighbouring element's son has already created a node at
1468 /// a position corresponding to the local fractional position within the
1469 /// present element, s_fraction, return
1470 /// a pointer to that node. If not, return NULL (0). If the node is
1471 /// periodic the flag is_periodic will be true
1472 //===================================================================
1473 template<unsigned INITIAL_NNODE_1D>
1476  bool &is_periodic)
1477 {
1478  using namespace QuadTreeNames;
1479 
1480  //Calculate the edges on which the node lies
1481  Vector<int> edges;
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);}
1486 
1487  //Find the number of edges
1488  unsigned n_size = edges.size();
1489  //If there are no edges, then there is no neighbour, return 0
1490  if(n_size==0) {return 0;}
1491 
1492  Vector<unsigned> translate_s(2);
1493  Vector<double> s_lo_neigh(2);
1494  Vector<double> s_hi_neigh(2);
1495  Vector<double> s(2);
1496 
1497  int neigh_edge, diff_level;
1498  bool in_neighbouring_tree;
1499 
1500  //Loop over the edges
1501  for(unsigned j=0;j<n_size;j++)
1502  {
1503  // Find pointer to neighbouring element along edge
1504  QuadTree* neigh_pt;
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);
1508 
1509  // Neighbour exists
1510  if(neigh_pt!=0)
1511  {
1512  // Have its nodes been created yet?
1513  // (Must look in sons of unfinished neighbours too!!!)
1514  if(1)
1515  {
1516  //We now need to translate the nodal location
1517  //as defined in terms of the fractional coordinates of the present
1518  //element into those of its neighbour
1519 
1520  //Calculate the local coordinate in the neighbour
1521  //Note that we need to use the translation scheme in case
1522  //the local coordinates are swapped in the neighbour.
1523  for(unsigned i=0;i<2;i++)
1524  {
1525  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
1526  (s_hi_neigh[i] - s_lo_neigh[i]);
1527  }
1528 
1529  // Check if the element has sons
1530  if(neigh_pt->nsons()!=0)
1531  {
1532  //First, find the son element in which the node should live
1533 
1534  //Find coordinates in the sons
1535  Vector<double> s_in_son(2);
1536  int son=-10;
1537  //If negative on the west side
1538  if(s[0] < 0.0)
1539  {
1540  //On the south side
1541  if(s[1] < 0.0)
1542  {
1543  //It's the southwest son
1544  son = SW;
1545  s_in_son[0] = 1.0 + 2.0*s[0];
1546  s_in_son[1] = 1.0 + 2.0*s[1];
1547  }
1548  //On the north side
1549  else
1550  {
1551  //It's the northwest son
1552  son = NW;
1553  s_in_son[0] = 1.0 + 2.0*s[0];
1554  s_in_son[1] = -1.0 + 2.0*s[1];
1555  }
1556  }
1557  else
1558  {
1559  //On the south side
1560  if(s[1] < 0.0)
1561  {
1562  //It's the southeast son
1563  son = SE;
1564  s_in_son[0] = -1.0 + 2.0*s[0];
1565  s_in_son[1] = 1.0 + 2.0*s[1];
1566  }
1567  //On the north side
1568  else
1569  {
1570  //It's the northeast son
1571  son = NE;
1572  s_in_son[0] = -1.0 + 2.0*s[0];
1573  s_in_son[1] = -1.0 + 2.0*s[1];
1574  }
1575  }
1576 
1577  //Find the node in the neighbour's son
1578  Node* neighbour_son_node_pt =
1579  neigh_pt->son_pt(son)->object_pt()->
1580  get_node_at_local_coordinate(s_in_son);
1581 
1582  //If there is a node, return it
1583  if(neighbour_son_node_pt!=0)
1584  {
1585  //Now work out whether it's a periodic boundary
1586  //only possible if we have moved into a neighbouring tree
1587  if(in_neighbouring_tree)
1588  {
1589  //Return whether the neighbour is periodic
1590  is_periodic =
1591  quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1592  }
1593  //Return the pointer to the neighbouring node
1594  return neighbour_son_node_pt;
1595  }
1596  }
1597  else
1598  {
1599  // No sons to search in, so no node can be found
1600  return 0;
1601  }
1602  }
1603  }
1604  }
1605  //Node not found, return null
1606  return 0;
1607 }
1608 
1609 //==================================================================
1610 /// Set the correct p-order of the element based on that of its
1611 /// father. Then construct an integration scheme of the correct order.
1612 /// If an adopted father is specified, information from this is
1613 /// used instead of using the father found from the tree.
1614 //==================================================================
1615 template<unsigned INITIAL_NNODE_1D>
1617 initial_setup(Tree* const &adopted_father_pt, const unsigned &initial_p_order)
1618 {
1619  //Storage for pointer to my father (in quadtree impersonation)
1620  QuadTree* father_pt = 0;
1621 
1622  // Check if an adopted father has been specified
1623  if (adopted_father_pt!=0)
1624  {
1625  //Get pointer to my father (in quadtree impersonation)
1626  father_pt = dynamic_cast<QuadTree*>(adopted_father_pt);
1627  }
1628  // Check if element is in a tree
1629  else if (Tree_pt!=0)
1630  {
1631  //Get pointer to my father (in quadtree impersonation)
1632  father_pt = dynamic_cast<QuadTree*>(quadtree_pt()->father_pt());
1633  }
1634  //else
1635  // {
1636  // throw OomphLibError(
1637  // "Element not in a tree, and no adopted father has been specified!",
1638  // OOMPH_CURRENT_FUNCTION,
1639  // OOMPH_EXCEPTION_LOCATION);
1640  // }
1641 
1642  // Check if we found a father
1643  if (father_pt!=0 || initial_p_order!=0)
1644  {
1645  if(father_pt!=0)
1646  {
1649  (father_pt->object_pt());
1650  if (father_el_pt!=0)
1651  {
1652  unsigned father_p_order = father_el_pt->p_order();
1653  // Set p-order to that of father
1654  P_order = father_p_order;
1655  }
1656  }
1657  else
1658  {
1659  P_order = initial_p_order;
1660  }
1661 
1662  // Now sort out the element...
1663  // (has p^2 nodes)
1664  unsigned new_n_node = P_order*P_order;
1665 
1666  // Allocate new space for Nodes (at the element level)
1667  this->set_n_node(new_n_node);
1668 
1669  // Set integration scheme
1670  delete this->integral_pt();
1671  switch(P_order)
1672  {
1673  case 2:
1674  this->set_integration_scheme(new GaussLobattoLegendre<2,2>);
1675  break;
1676  case 3:
1677  this->set_integration_scheme(new GaussLobattoLegendre<2,3>);
1678  break;
1679  case 4:
1680  this->set_integration_scheme(new GaussLobattoLegendre<2,4>);
1681  break;
1682  case 5:
1683  this->set_integration_scheme(new GaussLobattoLegendre<2,5>);
1684  break;
1685  case 6:
1686  this->set_integration_scheme(new GaussLobattoLegendre<2,6>);
1687  break;
1688  case 7:
1689  this->set_integration_scheme(new GaussLobattoLegendre<2,7>);
1690  break;
1691  default:
1692  std::ostringstream error_message;
1693  error_message <<"\nERROR: Exceeded maximum polynomial order for";
1694  error_message <<"\n integration scheme.\n";
1695  throw OomphLibError(error_message.str(),
1696  OOMPH_CURRENT_FUNCTION,
1697  OOMPH_EXCEPTION_LOCATION);
1698  }
1699  }
1700 }
1701 
1702 //==================================================================
1703 /// Check the father element for any required nodes which
1704 /// already exist
1705 //==================================================================
1706 template<unsigned INITIAL_NNODE_1D>
1708  Mesh*& mesh_pt,
1709  Vector<Node*>& new_node_pt)
1710 {
1711  using namespace QuadTreeNames;
1712 
1713  //Get the number of 1d nodes
1714  unsigned n_p = nnode_1d();
1715 
1716  //Check whether static father_bound needs to be created
1717  if(Father_bound[n_p].nrow()==0) {setup_father_bounds();}
1718 
1719  //Pointer to my father (in quadtree impersonation)
1720  QuadTree* father_pt = dynamic_cast<QuadTree*>(quadtree_pt()->father_pt());
1721 
1722  // What type of son am I? Ask my quadtree representation...
1723  int son_type = Tree_pt->son_type();
1724 
1725  // Has somebody build me already? (If any nodes have not been built)
1726  if (!nodes_built())
1727  {
1728 #ifdef PARANOID
1729  if (father_pt==0)
1730  {
1731  std::string error_message =
1732  "Something fishy here: I have no father and yet \n";
1733  error_message +=
1734  "I have no nodes. Who has created me then?!\n";
1735 
1736  throw OomphLibError(error_message,
1737  OOMPH_CURRENT_FUNCTION,
1738  OOMPH_EXCEPTION_LOCATION);
1739  }
1740 #endif
1741 
1742  // Return pointer to father element
1743  RefineableQElement<2>* father_el_pt
1744  = dynamic_cast<RefineableQElement<2>*>(father_pt->object_pt());
1745 
1746  // Timestepper should be the same for all nodes in father
1747  // element -- use it create timesteppers for new nodes
1748  TimeStepper* time_stepper_pt=father_el_pt->node_pt(0)->time_stepper_pt();
1749 
1750  // Number of history values (incl. present)
1751  unsigned ntstorage=time_stepper_pt->ntstorage();
1752 
1753  Vector<double> s_lo(2);
1754  Vector<double> s_hi(2);
1755  Vector<double> s(2);
1756  Vector<double> x(2);
1757 
1758  // Setup vertex coordinates in father element:
1759  //--------------------------------------------
1760  switch(son_type)
1761  {
1762  case SW:
1763  s_lo[0]=-1.0;
1764  s_hi[0]= 0.0;
1765  s_lo[1]=-1.0;
1766  s_hi[1]= 0.0;
1767  break;
1768 
1769  case SE:
1770  s_lo[0]= 0.0;
1771  s_hi[0]= 1.0;
1772  s_lo[1]=-1.0;
1773  s_hi[1]= 0.0;
1774  break;
1775 
1776  case NE:
1777  s_lo[0]= 0.0;
1778  s_hi[0]= 1.0;
1779  s_lo[1]= 0.0;
1780  s_hi[1]= 1.0;
1781  break;
1782 
1783  case NW:
1784  s_lo[0]=-1.0;
1785  s_hi[0]= 0.0;
1786  s_lo[1]= 0.0;
1787  s_hi[1]= 1.0;
1788  break;
1789  }
1790 
1791  // Pass macro element pointer on to sons and
1792  // set coordinates in macro element
1793  // hierher why can I see this?
1794  if(father_el_pt->macro_elem_pt()!=0)
1795  {
1796  set_macro_elem_pt(father_el_pt->macro_elem_pt());
1797  for(unsigned i=0;i<2;i++)
1798  {
1799  s_macro_ll(i)= father_el_pt->s_macro_ll(i)+
1800  0.5*(s_lo[i]+1.0)*(father_el_pt->s_macro_ur(i)-
1801  father_el_pt->s_macro_ll(i));
1802  s_macro_ur(i)= father_el_pt->s_macro_ll(i)+
1803  0.5*(s_hi[i]+1.0)*(father_el_pt->s_macro_ur(i)-
1804  father_el_pt->s_macro_ll(i));
1805  }
1806  }
1807 
1808 
1809  // If the father element hasn't been generated yet, we're stuck...
1810  if(father_el_pt->node_pt(0)==0)
1811  {
1812  throw OomphLibError(
1813  "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1814  OOMPH_CURRENT_FUNCTION,
1815  OOMPH_EXCEPTION_LOCATION);
1816  }
1817  else
1818  {
1819  unsigned jnod=0;
1820  Vector<double> x_small(2);
1821  Vector<double> x_large(2);
1822 
1823  Vector<double> s_fraction(2);
1824  // Loop over nodes in element
1825  for(unsigned i0=0;i0<n_p;i0++)
1826  {
1827  //Get the fractional position of the node in the direction of s[0]
1828  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
1829  // Local coordinate in father element
1830  s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
1831 
1832  for(unsigned i1=0;i1<n_p;i1++)
1833  {
1834  //Get the fractional position of the node in the direction of s[1]
1835  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
1836  // Local coordinate in father element
1837  s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
1838 
1839  // Local node number
1840  jnod= i0 + n_p*i1;
1841 
1842  //Check whether the father's node is periodic if so, complain
1843  /* {
1844  Node* father_node_pt = father_el_pt->node_pt(jnod);
1845  if((father_node_pt->is_a_copy()) ||
1846  (father_node_pt->position_is_a_copy()))
1847  {
1848  throw OomphLibError(
1849  "Can't handle periodic nodes (yet).",
1850  OOMPH_CURRENT_FUNCTION,
1851  OOMPH_EXCEPTION_LOCATION);
1852  }
1853  }*/
1854 
1855  // Initialise flag: So far, this node hasn't been built
1856  // or copied yet
1857  //bool node_done=false;
1858 
1859  //Get the pointer to the node in the father, returns NULL
1860  //if there is not node
1861  Node* created_node_pt = father_el_pt->get_node_at_local_coordinate(s);
1862 
1863  // Does this node already exist in father element?
1864  //------------------------------------------------
1865  if(created_node_pt!=0)
1866  {
1867  // Copy node across
1868  this->node_pt(jnod) = created_node_pt;
1869 
1870  //Make sure that we update the values at the node so that
1871  //they are consistent with the present representation.
1872  //This is only need for mixed interpolation where the value
1873  //at the father could now become active.
1874 
1875  // Loop over all history values
1876  for(unsigned t=0;t<ntstorage;t++)
1877  {
1878  // Get values from father element
1879  // Note: get_interpolated_values() sets Vector size itself.
1880  Vector<double> prev_values;
1881  father_el_pt->get_interpolated_values(t,s,prev_values);
1882  //Find the minimum number of values
1883  //(either those stored at the node, or those returned by
1884  // the function)
1885  unsigned n_val_at_node = created_node_pt->nvalue();
1886  unsigned n_val_from_function = prev_values.size();
1887  //Use the ternary conditional operator here
1888  unsigned n_var = n_val_at_node < n_val_from_function ?
1889  n_val_at_node : n_val_from_function;
1890  //Assign the values that we can
1891  for(unsigned k=0;k<n_var;k++)
1892  {
1893  created_node_pt->set_value(t,k,prev_values[k]);
1894  }
1895  }
1896 
1897  // Node has been created by copy
1898  //node_done=true;
1899  }
1900  } // End of vertical loop over nodes in element
1901  } // End of horizontal loop over nodes in element
1902  } // Sanity check: Father element has been generated
1903 
1904  } // End of nodes not built
1905  else
1906  {
1907  // If this is element updates by macro element node update then we need
1908  // to set the correct info in the nodes here since this won't get done
1909  // later in build() because we already have all our nodes created.
1910  MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
1912  if (m_el_pt!=0)
1913  {
1914  // Get vector of geometric objects
1915  Vector<GeomObject*> geom_object_pt = m_el_pt->geom_object_pt();
1916 
1917  //// Build update info by passing vector of geometric objects:
1918  //// This sets the current element to be the update element
1919  //// for all of the element's nodes -- this is reversed
1920  //// if the element is ever un-refined in the father element's
1921  //// rebuild_from_sons() function which overwrites this
1922  //// assignment to avoid nasty segmentation faults that occur
1923  //// when a node tries to update itself via an element that no
1924  //// longer exists...
1925  m_el_pt->set_node_update_info(geom_object_pt);
1926  }
1927  }
1928 }
1929 
1930 //==================================================================
1931 /// p-refine the element inc times. (If inc<0 then p-unrefinement
1932 /// is performed.)
1933 //==================================================================
1934 template<unsigned INITIAL_NNODE_1D>
1936  const int &inc,
1937  Mesh* const &mesh_pt,
1938  GeneralisedElement* const &clone_pt)
1939 {
1940  //Cast clone to correct type
1942  = dynamic_cast<PRefineableQElement<2,INITIAL_NNODE_1D>*>(clone_pt);
1943 
1944  // If it is a MacroElement then we need to copy the geometric data too.
1946  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
1947  if(m_el_pt!=0)
1948  {
1949  MacroElementNodeUpdateElementBase* clone_m_el_pt =
1950  dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_pt);
1951 
1952  // Store local copy of geom object vector, so it can be passed on
1953  // to son elements (and their nodes) during refinement
1954  unsigned ngeom_object=m_el_pt->ngeom_object();
1955  clone_m_el_pt->geom_object_pt().resize(ngeom_object);
1956  for (unsigned i=0;i<ngeom_object;i++)
1957  {
1958  clone_m_el_pt->geom_object_pt()[i]=m_el_pt->geom_object_pt(i);
1959  }
1960 
1961  //clone_m_el_pt->set_node_update_info(m_el_pt->geom_object_pt());
1962  }
1963 
1964  //Check if we can use it
1965  if(clone_el_pt==0)
1966  {
1967  throw OomphLibError(
1968  "Cloned copy must be a PRefineableQElement<2,INITIAL_NNODE_1D>!",
1969  OOMPH_CURRENT_FUNCTION,
1970  OOMPH_EXCEPTION_LOCATION);
1971  }
1972 #ifdef PARANOID
1973  //Clone exists, so check that it is infact a copy of me
1974  else
1975  {
1976  //Flag to keep track of check
1977  bool clone_is_ok = true;
1978 
1979  //Does the clone have the correct p-order?
1980  clone_is_ok = clone_is_ok && (clone_el_pt->p_order() == this->p_order());
1981 
1982  if(!clone_is_ok)
1983  {
1984  std::ostringstream err_stream;
1985  err_stream << "Clone element has a different p-order from me,"
1986  << std::endl
1987  << "but it is supposed to be a copy!" << std::endl;
1988 
1989  throw OomphLibError(err_stream.str(),
1990  OOMPH_CURRENT_FUNCTION,
1991  OOMPH_EXCEPTION_LOCATION);
1992  }
1993 
1994  //Does the clone have the same number of nodes as me?
1995  clone_is_ok = clone_is_ok && (clone_el_pt->nnode() == this->nnode());
1996 
1997  if(!clone_is_ok)
1998  {
1999  std::ostringstream err_stream;
2000  err_stream << "Clone element has a different number of nodes from me,"
2001  << std::endl
2002  << "but it is supposed to be a copy!" << std::endl;
2003 
2004  throw OomphLibError(err_stream.str(),
2005  OOMPH_CURRENT_FUNCTION,
2006  OOMPH_EXCEPTION_LOCATION);
2007  }
2008 
2009  //Does the clone have the same nodes as me?
2010  for(unsigned n=0; n<this->nnode(); n++)
2011  {
2012  clone_is_ok = clone_is_ok && (clone_el_pt->node_pt(n) == this->node_pt(n));
2013  }
2014 
2015  if(!clone_is_ok)
2016  {
2017  std::ostringstream err_stream;
2018  err_stream << "The nodes belonging to the clone element are different"
2019  << std::endl
2020  << "from mine, but it is supposed to be a copy!" << std::endl;
2021 
2022  throw OomphLibError(err_stream.str(),
2023  OOMPH_CURRENT_FUNCTION,
2024  OOMPH_EXCEPTION_LOCATION);
2025  }
2026 
2027  // Is this a macro element?
2029  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
2030  if(m_el_pt!=0)
2031  {
2032  // Get macro element version of clone
2033  MacroElementNodeUpdateElementBase* clone_m_el_pt =
2034  dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_el_pt);
2035 
2036  //Does the clone have the same geometric objects?
2037  Vector<GeomObject*> clone_geom_object_pt(clone_m_el_pt->geom_object_pt());
2038  Vector<GeomObject*> geom_object_pt(m_el_pt->geom_object_pt());
2039 
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++)
2042  {
2043  clone_is_ok = clone_is_ok && (clone_geom_object_pt[n] == geom_object_pt[n]);
2044  }
2045 
2046  if(!clone_is_ok)
2047  {
2048  std::ostringstream err_stream;
2049  err_stream << "The clone element has different geometric objects"
2050  << std::endl
2051  << "from mine, but it is supposed to be a copy!" << std::endl;
2052 
2053  throw OomphLibError(err_stream.str(),
2054  "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2055  OOMPH_EXCEPTION_LOCATION);
2056  }
2057  }
2058 
2059  //If we get to here then the clone has all the information we require
2060  }
2061 #endif
2062 
2063  // Currently we can't handle the case of generalised coordinates
2064  // since we haven't established how they should be interpolated.
2065  // Buffer this case:
2066  if (clone_el_pt->node_pt(0)->nposition_type()!=1)
2067  {
2068  throw OomphLibError(
2069  "Can't handle generalised nodal positions (yet).",
2070  OOMPH_CURRENT_FUNCTION,
2071  OOMPH_EXCEPTION_LOCATION);
2072  }
2073 
2074  // Timestepper should be the same for all nodes -- use it
2075  // to create timesteppers for new nodes
2076  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
2077 
2078  // Get number of history values (incl. present)
2079  unsigned ntstorage = time_stepper_pt->ntstorage();
2080 
2081  // Increment p-order of the element
2082  P_order += inc;
2083 
2084  // Change integration scheme
2085  delete this->integral_pt();
2086  switch(P_order)
2087  {
2088  case 2:
2089  this->set_integration_scheme(new GaussLobattoLegendre<2,2>);
2090  break;
2091  case 3:
2092  this->set_integration_scheme(new GaussLobattoLegendre<2,3>);
2093  break;
2094  case 4:
2095  this->set_integration_scheme(new GaussLobattoLegendre<2,4>);
2096  break;
2097  case 5:
2098  this->set_integration_scheme(new GaussLobattoLegendre<2,5>);
2099  break;
2100  case 6:
2101  this->set_integration_scheme(new GaussLobattoLegendre<2,6>);
2102  break;
2103  case 7:
2104  this->set_integration_scheme(new GaussLobattoLegendre<2,7>);
2105  break;
2106  default:
2107  std::ostringstream error_message;
2108  error_message <<"\nERROR: Exceeded maximum polynomial order for";
2109  error_message <<"\n integration scheme.\n";
2110  throw OomphLibError(error_message.str(),
2111  OOMPH_CURRENT_FUNCTION,
2112  OOMPH_EXCEPTION_LOCATION);
2113  }
2114 
2115  // Allocate new space for Nodes (at the element level)
2116  this->set_n_node(P_order*P_order);
2117 
2118  // Copy vertex nodes and create new edge and internal nodes
2119  //---------------------------------------------------------
2120 
2121  // Setup vertex coordinates in element:
2122  //-------------------------------------
2123  Vector<double> s_lo(2);
2124  Vector<double> s_hi(2);
2125  s_lo[0]=-1.0;
2126  s_hi[0]= 1.0;
2127  s_lo[1]=-1.0;
2128  s_hi[1]= 1.0;
2129 
2130  //Local coordinate in element
2131  Vector<double> s(2);
2132 
2133  unsigned jnod=0;
2134 
2135  Vector<double> s_fraction(2);
2136  // Loop over nodes in element
2137  for(unsigned i0=0;i0<P_order;i0++)
2138  {
2139  //Get the fractional position of the node in the direction of s[0]
2140  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
2141  // Local coordinate
2142  s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
2143 
2144  for(unsigned i1=0;i1<P_order;i1++)
2145  {
2146  //Get the fractional position of the node in the direction of s[1]
2147  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
2148  // Local coordinate
2149  s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
2150 
2151  // Local node number
2152  jnod= i0 + P_order*i1;
2153 
2154  // Initialise flag: So far, this node hasn't been built
2155  // or copied yet
2156  bool node_done=false;
2157 
2158  //Get the pointer to the node in this element (or rather, its clone),
2159  //returns NULL if there is not node
2160  Node* created_node_pt = clone_el_pt->get_node_at_local_coordinate(s);
2161 
2162  // Does this node already exist in this element?
2163  //----------------------------------------------
2164  if (created_node_pt!=0)
2165  {
2166  // Copy node across
2167  this->node_pt(jnod) = created_node_pt;
2168 
2169  //Make sure that we update the values at the node so that
2170  //they are consistent with the present representation.
2171  //This is only need for mixed interpolation where the value
2172  //at the father could now become active.
2173 
2174  // Loop over all history values
2175  for(unsigned t=0;t<ntstorage;t++)
2176  {
2177  // Get values from father element
2178  // Note: get_interpolated_values() sets Vector size itself.
2179  Vector<double> prev_values;
2180  clone_el_pt->get_interpolated_values(t,s,prev_values);
2181  //Find the minimum number of values
2182  //(either those stored at the node, or those returned by
2183  // the function)
2184  unsigned n_val_at_node = created_node_pt->nvalue();
2185  unsigned n_val_from_function = prev_values.size();
2186  //Use the ternary conditional operator here
2187  unsigned n_var = n_val_at_node < n_val_from_function ?
2188  n_val_at_node : n_val_from_function;
2189  //Assign the values that we can
2190  for(unsigned k=0;k<n_var;k++)
2191  {
2192  created_node_pt->set_value(t,k,prev_values[k]);
2193  }
2194  }
2195 
2196  // Node has been created by copy
2197  node_done = true;
2198  }
2199  // Node does not exist in this element but might already
2200  //------------------------------------------------------
2201  // have been created by neighbouring elements
2202  //-------------------------------------------
2203  else
2204  {
2205  //Was the node created by one of its neighbours
2206  //Whether or not the node lies on an edge can be calculated
2207  //by from the fractional position
2208  bool is_periodic = false;
2209  created_node_pt =
2210  node_created_by_neighbour(s_fraction,is_periodic);
2211 
2212  //If the node was so created, assign the pointers
2213  if (created_node_pt!=0)
2214  {
2215  //If the node is periodic
2216  if(is_periodic)
2217  {
2218  //Now the node must be on a boundary, but we don't know which
2219  //one
2220  //The returned created_node_pt is actually the neighbouring
2221  //periodic node
2222  Node* neighbour_node_pt = created_node_pt;
2223 
2224  // Determine the edge on which the new node will live
2225  //(cannot be a vertex node)
2226  int my_bound = Tree::OMEGA;
2227  if(s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2228  else if(s_fraction[0] == 1.0) my_bound = QuadTreeNames::E;
2229  else if(s_fraction[1] == 0.0) my_bound = QuadTreeNames::S;
2230  else if(s_fraction[1] == 1.0) my_bound = QuadTreeNames::N;
2231 
2232  // Storage for the set of Mesh boundaries on which the
2233  // appropriate edge lives.
2234  // [New nodes should always be mid-edge nodes and therefore
2235  //only live on one boundary but just to play it safe...]
2236  std::set<unsigned> boundaries;
2237  //Only get the boundaries if we are at the edge of
2238  //an element. Nodes in the centre of an element cannot be
2239  //on Mesh boundaries
2240  if(my_bound!=Tree::OMEGA)
2241  {clone_el_pt->get_boundaries(my_bound,boundaries);}
2242 
2243 #ifdef PARANOID
2244  //Case where a new node lives on more than one boundary
2245  // seems fishy enough to flag
2246  if (boundaries.size()>1)
2247  {
2248  throw OomphLibError(
2249  "boundaries.size()!=1 seems a bit strange..\n",
2250  OOMPH_CURRENT_FUNCTION,
2251  OOMPH_EXCEPTION_LOCATION);
2252  }
2253 
2254  //Case when there are no boundaries, we are in big trouble
2255  if(boundaries.size() == 0)
2256  {
2257  std::ostringstream error_stream;
2258  error_stream
2259  << "Periodic node is not on a boundary...\n"
2260  << "Coordinates: "
2261  << created_node_pt->x(0) << " "
2262  << created_node_pt->x(1) << "\n";
2263  throw OomphLibError(
2264  error_stream.str(),
2265  OOMPH_CURRENT_FUNCTION,
2266  OOMPH_EXCEPTION_LOCATION);
2267  }
2268 #endif
2269 
2270  // Create node and set the pointer to it from the element
2271  created_node_pt =
2272  this->construct_boundary_node(jnod,time_stepper_pt);
2273  //Make the node periodic from the neighbour
2274  created_node_pt->
2275  make_periodic(neighbour_node_pt);
2276 
2277  // Loop over # of history values
2278  for (unsigned t=0;t<ntstorage;t++)
2279  {
2280  // Get position from father element -- this uses the macro
2281  // element representation if appropriate. If the node
2282  // turns out to be a hanging node later on, then
2283  // its position gets adjusted in line with its
2284  // hanging node interpolation.
2285  Vector<double> x_prev(2);
2286  clone_el_pt->get_x(t,s,x_prev);
2287  // Set previous positions of the new node
2288  for(unsigned i=0;i<2;i++)
2289  {
2290  created_node_pt->x(t,i) = x_prev[i];
2291  }
2292  }
2293 
2294  // Check if we need to add nodes to the mesh
2295  if(mesh_pt!=0)
2296  {
2297  // Next, we Update the boundary lookup schemes
2298  //Loop over the boundaries stored in the set
2299  for(std::set<unsigned>::iterator it = boundaries.begin();
2300  it != boundaries.end(); ++it)
2301  {
2302  //Add the node to the boundary
2303  mesh_pt->add_boundary_node(*it,created_node_pt);
2304 
2305  //If we have set an intrinsic coordinate on this
2306  //mesh boundary then it must also be interpolated on
2307  //the new node
2308  //Now interpolate the intrinsic boundary coordinate
2309  if(mesh_pt->boundary_coordinate_exists(*it)==true)
2310  {
2311  Vector<double> zeta(1);
2312  clone_el_pt->interpolated_zeta_on_edge(*it,
2313  my_bound,
2314  s,zeta);
2315 
2316  created_node_pt->set_coordinates_on_boundary(*it,zeta);
2317  }
2318  }
2319 
2320  //Make sure that we add the node to the mesh
2321  mesh_pt->add_node_pt(created_node_pt);
2322  }
2323  } //End of periodic case
2324  //Otherwise the node is not periodic, so just set the
2325  //pointer to the neighbours node
2326  else
2327  {
2328  this->node_pt(jnod) = created_node_pt;
2329  }
2330  //Node has been created
2331  node_done = true;
2332  }
2333  // Node does not exist in neighbour element but might already
2334  //-----------------------------------------------------------
2335  // have been created by a son of a neighbouring element
2336  //-----------------------------------------------------
2337  else
2338  {
2339  //Was the node created by one of its neighbours' sons
2340  //Whether or not the node lies on an edge can be calculated
2341  //by from the fractional position
2342  bool is_periodic = false;
2343  created_node_pt =
2344  node_created_by_son_of_neighbour(s_fraction,is_periodic);
2345 
2346  //If the node was so created, assign the pointers
2347  if (created_node_pt!=0)
2348  {
2349  //If the node is periodic
2350  if(is_periodic)
2351  {
2352  //Now the node must be on a boundary, but we don't know which
2353  //one
2354  //The returned created_node_pt is actually the neighbouring
2355  //periodic node
2356  Node* neighbour_node_pt = created_node_pt;
2357 
2358  // Determine the edge on which the new node will live
2359  //(cannot be a vertex node)
2360  int my_bound = Tree::OMEGA;
2361  if(s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2362  else if(s_fraction[0] == 1.0) my_bound = QuadTreeNames::E;
2363  else if(s_fraction[1] == 0.0) my_bound = QuadTreeNames::S;
2364  else if(s_fraction[1] == 1.0) my_bound = QuadTreeNames::N;
2365 
2366  // Storage for the set of Mesh boundaries on which the
2367  // appropriate edge lives.
2368  // [New nodes should always be mid-edge nodes and therefore
2369  //only live on one boundary but just to play it safe...]
2370  std::set<unsigned> boundaries;
2371  //Only get the boundaries if we are at the edge of
2372  //an element. Nodes in the centre of an element cannot be
2373  //on Mesh boundaries
2374  if(my_bound!=Tree::OMEGA)
2375  {clone_el_pt->get_boundaries(my_bound,boundaries);}
2376 
2377 #ifdef PARANOID
2378  //Case where a new node lives on more than one boundary
2379  // seems fishy enough to flag
2380  if (boundaries.size()>1)
2381  {
2382  throw OomphLibError(
2383  "boundaries.size()!=1 seems a bit strange..\n",
2384  OOMPH_CURRENT_FUNCTION,
2385  OOMPH_EXCEPTION_LOCATION);
2386  }
2387 
2388  //Case when there are no boundaries, we are in big trouble
2389  if(boundaries.size() == 0)
2390  {
2391  std::ostringstream error_stream;
2392  error_stream
2393  << "Periodic node is not on a boundary...\n"
2394  << "Coordinates: "
2395  << created_node_pt->x(0) << " "
2396  << created_node_pt->x(1) << "\n";
2397  throw OomphLibError(
2398  error_stream.str(),
2399  OOMPH_CURRENT_FUNCTION,
2400  OOMPH_EXCEPTION_LOCATION);
2401  }
2402 #endif
2403 
2404  // Create node and set the pointer to it from the element
2405  created_node_pt =
2406  this->construct_boundary_node(jnod,time_stepper_pt);
2407  //Make the node periodic from the neighbour
2408  created_node_pt->
2409  make_periodic(neighbour_node_pt);
2410 
2411  // Loop over # of history values
2412  for (unsigned t=0;t<ntstorage;t++)
2413  {
2414  // Get position from father element -- this uses the macro
2415  // element representation if appropriate. If the node
2416  // turns out to be a hanging node later on, then
2417  // its position gets adjusted in line with its
2418  // hanging node interpolation.
2419  Vector<double> x_prev(2);
2420  clone_el_pt->get_x(t,s,x_prev);
2421  // Set previous positions of the new node
2422  for(unsigned i=0;i<2;i++)
2423  {
2424  created_node_pt->x(t,i) = x_prev[i];
2425  }
2426  }
2427 
2428  // Check if we need to add nodes to the mesh
2429  if(mesh_pt!=0)
2430  {
2431  // Next, we Update the boundary lookup schemes
2432  //Loop over the boundaries stored in the set
2433  for(std::set<unsigned>::iterator it = boundaries.begin();
2434  it != boundaries.end(); ++it)
2435  {
2436  //Add the node to the boundary
2437  mesh_pt->add_boundary_node(*it,created_node_pt);
2438 
2439  //If we have set an intrinsic coordinate on this
2440  //mesh boundary then it must also be interpolated on
2441  //the new node
2442  //Now interpolate the intrinsic boundary coordinate
2443  if(mesh_pt->boundary_coordinate_exists(*it)==true)
2444  {
2445  Vector<double> zeta(1);
2446  clone_el_pt->interpolated_zeta_on_edge(*it,
2447  my_bound,
2448  s,zeta);
2449 
2450  created_node_pt->set_coordinates_on_boundary(*it,zeta);
2451  }
2452  }
2453 
2454  //Make sure that we add the node to the mesh
2455  mesh_pt->add_node_pt(created_node_pt);
2456  }
2457  } //End of periodic case
2458  //Otherwise the node is not periodic, so just set the
2459  //pointer to the neighbours node
2460  else
2461  {
2462  this->node_pt(jnod) = created_node_pt;
2463  }
2464  //Node has been created
2465  node_done = true;
2466  } //Node does not exist in son of neighbouring element
2467  } //Node does not exist in neighbouring element
2468  } // Node does not exist in this element
2469 
2470  // Node has not been built anywhere ---> build it here
2471  if (!node_done)
2472  {
2473  //Firstly, we need to determine whether or not a node lies
2474  //on the boundary before building it, because
2475  //we actually assign a different type of node on boundaries.
2476 
2477  // Determine the edge on which the new node will live
2478  //(cannot be a vertex node)
2479  int my_bound = Tree::OMEGA;
2480  if(s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2481  else if(s_fraction[0] == 1.0) my_bound = QuadTreeNames::E;
2482  else if(s_fraction[1] == 0.0) my_bound = QuadTreeNames::S;
2483  else if(s_fraction[1] == 1.0) my_bound = QuadTreeNames::N;
2484 
2485  // Storage for the set of Mesh boundaries on which the
2486  // appropriate edge lives.
2487  // [New nodes should always be mid-edge nodes and therefore
2488  //only live on one boundary but just to play it safe...]
2489  std::set<unsigned> boundaries;
2490  //Only get the boundaries if we are at the edge of
2491  //an element. Nodes in the centre of an element cannot be
2492  //on Mesh boundaries
2493  if(my_bound!=Tree::OMEGA)
2494  {clone_el_pt->get_boundaries(my_bound,boundaries);}
2495 
2496 #ifdef PARANOID
2497  //Case where a new node lives on more than one boundary
2498  // seems fishy enough to flag
2499  if (boundaries.size()>1)
2500  {
2501  throw OomphLibError(
2502  "boundaries.size()!=1 seems a bit strange..\n",
2503  OOMPH_CURRENT_FUNCTION,
2504  OOMPH_EXCEPTION_LOCATION);
2505  }
2506 #endif
2507 
2508  //If the node lives on a mesh boundary,
2509  //then we need to create a boundary node
2510  if(boundaries.size()> 0)
2511  {
2512  // Create node and set the pointer to it from the element
2513  created_node_pt = this->construct_boundary_node(jnod,time_stepper_pt);
2514 
2515  //Now we need to work out whether to pin the values at
2516  //the new node based on the boundary conditions applied at
2517  //its Mesh boundary
2518 
2519  //Get the boundary conditions from the father
2520  Vector<int> bound_cons(this->ncont_interpolated_values());
2521  clone_el_pt->get_bcs(my_bound,bound_cons);
2522 
2523  //Loop over the values and pin, if necessary
2524  unsigned n_value = created_node_pt->nvalue();
2525  for(unsigned k=0;k<n_value;k++)
2526  {
2527  if (bound_cons[k]) {created_node_pt->pin(k);}
2528  }
2529 
2530  // Solid node? If so, deal with the positional boundary
2531  // conditions:
2532  SolidNode* solid_node_pt =
2533  dynamic_cast<SolidNode*>(created_node_pt);
2534  if (solid_node_pt!=0)
2535  {
2536  //Get the positional boundary conditions from the father:
2537  unsigned n_dim = created_node_pt->ndim();
2538  Vector<int> solid_bound_cons(n_dim);
2539  RefineableSolidQElement<2>* clone_solid_el_pt=
2540  dynamic_cast<RefineableSolidQElement<2>*>(clone_el_pt);
2541 #ifdef PARANOID
2542  if (clone_solid_el_pt==0)
2543  {
2544  std::string error_message =
2545  "We have a SolidNode outside a refineable SolidElement\n";
2546  error_message +=
2547  "during mesh refinement -- this doesn't make sense";
2548 
2549  throw OomphLibError(
2550  error_message,
2551  OOMPH_CURRENT_FUNCTION,
2552  OOMPH_EXCEPTION_LOCATION);
2553  }
2554 #endif
2555  clone_solid_el_pt->
2556  get_solid_bcs(my_bound,solid_bound_cons);
2557 
2558  //Loop over the positions and pin, if necessary
2559  for(unsigned k=0;k<n_dim;k++)
2560  {
2561  if (solid_bound_cons[k]) {solid_node_pt->pin_position(k);}
2562  }
2563  } //End of if solid_node_pt
2564 
2565 
2566 
2567  // Check if we need to add nodes to the mesh
2568  if(mesh_pt!=0)
2569  {
2570  // Next, we Update the boundary lookup schemes
2571  //Loop over the boundaries stored in the set
2572  for(std::set<unsigned>::iterator it = boundaries.begin();
2573  it != boundaries.end(); ++it)
2574  {
2575  //Add the node to the boundary
2576  mesh_pt->add_boundary_node(*it,created_node_pt);
2577 
2578  //If we have set an intrinsic coordinate on this
2579  //mesh boundary then it must also be interpolated on
2580  //the new node
2581  //Now interpolate the intrinsic boundary coordinate
2582  if(mesh_pt->boundary_coordinate_exists(*it)==true)
2583  {
2584  Vector<double> zeta(1);
2585  clone_el_pt->interpolated_zeta_on_edge(*it,
2586  my_bound,
2587  s,zeta);
2588 
2589  created_node_pt->set_coordinates_on_boundary(*it,zeta);
2590  }
2591  }
2592  }
2593  }
2594  //Otherwise the node is not on a Mesh boundary and
2595  //we create a normal "bulk" node
2596  else
2597  {
2598  // Create node and set the pointer to it from the element
2599  created_node_pt = this->construct_node(jnod,time_stepper_pt);
2600  }
2601 
2602  //Now we set the position and values at the newly created node
2603 
2604  // In the first instance use macro element or FE representation
2605  // to create past and present nodal positions.
2606  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
2607  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
2608  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
2609  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
2610  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
2611  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
2612  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
2613 
2614  // Loop over # of history values
2615  for (unsigned t=0;t<ntstorage;t++)
2616  {
2617  // Get position from father element -- this uses the macro
2618  // element representation if appropriate. If the node
2619  // turns out to be a hanging node later on, then
2620  // its position gets adjusted in line with its
2621  // hanging node interpolation.
2622  Vector<double> x_prev(2);
2623  clone_el_pt->get_x(t,s,x_prev);
2624 
2625  // Set previous positions of the new node
2626  for(unsigned i=0;i<2;i++)
2627  {
2628  created_node_pt->x(t,i) = x_prev[i];
2629  }
2630  }
2631 
2632  // Loop over all history values
2633  for (unsigned t=0;t<ntstorage;t++)
2634  {
2635  // Get values from father element
2636  // Note: get_interpolated_values() sets Vector size itself.
2637  Vector<double> prev_values;
2638  clone_el_pt->get_interpolated_values(t,s,prev_values);
2639  //Initialise the values at the new node
2640  unsigned n_value = created_node_pt->nvalue();
2641  for(unsigned k=0;k<n_value;k++)
2642  {
2643  created_node_pt->set_value(t,k,prev_values[k]);
2644  }
2645  }
2646 
2647  // Add new node to mesh (if requested)
2648  if(mesh_pt!=0)
2649  {
2650  mesh_pt->add_node_pt(created_node_pt);
2651  }
2652 
2653  AlgebraicElementBase* alg_el_pt=
2654  dynamic_cast<AlgebraicElementBase*>(this);
2655 
2656  //If we do have an algebraic element
2657  if(alg_el_pt!=0)
2658  {
2659  std::string error_message =
2660  "Have not implemented p-refinement for";
2661  error_message +=
2662  "Algebraic p-refineable elements yet\n";
2663 
2664  throw
2665  OomphLibError(error_message,
2666  "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2667  OOMPH_EXCEPTION_LOCATION);
2668  }
2669 
2670  } //End of case when we build the node ourselves
2671 
2672  // Check if the element is an algebraic element
2673  AlgebraicElementBase* alg_el_pt=
2674  dynamic_cast<AlgebraicElementBase*>(this);
2675 
2676  // If the element is an algebraic element, setup
2677  // node position (past and present) from algebraic node update
2678  // function. This over-writes previous assingments that
2679  // were made based on the macro-element/FE representation.
2680  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
2681  // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
2682  // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
2683  // INFO FOR *ALL* ROOT ELEMENTS!
2684  if (alg_el_pt!=0)
2685  {
2686  // Build algebraic node update info for new node
2687  // This sets up the node update data for all node update
2688  // functions that are shared by all nodes in the father
2689  // element
2690  alg_el_pt->setup_algebraic_node_update(this->node_pt(jnod),s,
2691  clone_el_pt);
2692  }
2693 
2694  } // End of vertical loop over nodes in element
2695 
2696  } // End of horizontal loop over nodes in element
2697 
2698 
2699  // If the element is a MacroElementNodeUpdateElement, set
2700  // the update parameters for the current element's nodes --
2701  // all this needs is the vector of (pointers to the)
2702  // geometric objects that affect the MacroElement-based
2703  // node update -- this needs to be done to set the node
2704  // update info for newly created nodes
2705  MacroElementNodeUpdateElementBase* clone_m_el_pt=dynamic_cast<
2706  MacroElementNodeUpdateElementBase*>(clone_el_pt);
2707  if (clone_m_el_pt!=0)
2708  {
2709  // Get vector of geometric objects from father (construct vector
2710  // via copy operation)
2711  Vector<GeomObject*> geom_object_pt(clone_m_el_pt->geom_object_pt());
2712 
2713  // Cast current element to MacroElementNodeUpdateElement:
2714  MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
2716 
2717 #ifdef PARANOID
2718  if (m_el_pt==0)
2719  {
2720  std::string error_message =
2721  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
2722  error_message +=
2723  "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
2724  error_message += "then I should be too....\n";
2725 
2726  throw OomphLibError(error_message,
2727  OOMPH_CURRENT_FUNCTION,
2728  OOMPH_EXCEPTION_LOCATION);
2729  }
2730 #endif
2731  // Build update info by passing vector of geometric objects:
2732  // This sets the current element to be the update element
2733  // for all of the element's nodes -- this is reversed
2734  // if the element is ever un-refined in the father element's
2735  // rebuild_from_sons() function which overwrites this
2736  // assignment to avoid nasty segmentation faults that occur
2737  // when a node tries to update itself via an element that no
2738  // longer exists...
2739  m_el_pt->set_node_update_info(geom_object_pt);
2740 
2741  //// Now loop over nodes in element
2742  //unsigned n_node = this->nnode();
2743  //for (unsigned j=0;j<n_node;j++)
2744  // {
2745  // // Get local coordinate in element (Vector sets its own size)
2746  // Vector<double> s_in_node_update_element;
2747  // this->local_coordinate_of_node(j,s_in_node_update_element);
2748  //
2749  // // Pass the lot to the node
2750  // static_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))->
2751  // set_node_update_info(this,s_in_node_update_element,m_el_pt->geom_object_pt());
2752  // }
2753 
2754  ////BENFLAG:
2755  //std::cout << "Checking that all the nodes have this as their update element..." << std::endl;
2756  ////std::cout << "this = " << this << std::endl;
2757  //for(unsigned j=0; j<this->nnode(); j++)
2758  // {
2759  // //std::cout << this->node_pt(j) << ": [" << this->node_pt(j)->x(0) << "," << this->node_pt(j)->x(1) << "] update element: " << dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))->node_update_element_pt() << std::endl;
2760  // MacroElementNodeUpdateNode* mac_nod_pt = dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j));
2761  // if(mac_nod_pt->node_update_element_pt()!=this)
2762  // {
2763  // std::cout << "Something's not right! Update element is wrong..." << std::endl;
2764  // }
2765  // FiniteElement* up_el_pt = dynamic_cast<FiniteElement*>(mac_nod_pt->node_update_element_pt());
2766  // bool not_good = true;
2767  // for(unsigned l=0; l<up_el_pt->nnode(); l++)
2768  // {
2769  // if(up_el_pt->node_pt(l)==mac_nod_pt)
2770  // {
2771  // not_good = false;
2772  // break;
2773  // }
2774  // }
2775  // if(not_good==true)
2776  // {
2777  // std::cout << "Macro node doesn't belong to its update element!" << std::endl;
2778  // }
2779  // }
2780 
2781 
2782  // Loop over all nodes in element again, to re-set the positions
2783  // This must be done using the new element's macro-element
2784  // representation, rather than the old version which may be
2785  // of a different p-order!
2786  for(unsigned i0=0;i0<P_order;i0++)
2787  {
2788  //Get the fractional position of the node in the direction of s[0]
2789  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
2790  // Local coordinate
2791  s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
2792 
2793  for(unsigned i1=0;i1<P_order;i1++)
2794  {
2795  //Get the fractional position of the node in the direction of s[1]
2796  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
2797  // Local coordinate
2798  s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
2799 
2800  // Local node number
2801  jnod= i0 + P_order*i1;
2802 
2803  // Loop over # of history values
2804  for (unsigned t=0;t<ntstorage;t++)
2805  {
2806  // Get position from father element -- this uses the macro
2807  // element representation if appropriate. If the node
2808  // turns out to be a hanging node later on, then
2809  // its position gets adjusted in line with its
2810  // hanging node interpolation.
2811  Vector<double> x_prev(2);
2812  this->get_x(t,s,x_prev);
2813 
2814  // Set previous positions of the new node
2815  for(unsigned i=0;i<2;i++)
2816  {
2817  this->node_pt(jnod)->x(t,i) = x_prev[i];
2818  }
2819  }
2820  }
2821  }
2822  }
2823 
2824  // Not necessary to delete the old nodes since all original nodes are in the
2825  // current mesh and so will be pruned as part of the mesh adaption process.
2826 
2827  // Do any further-build required
2828  this->further_build();
2829 }
2830 
2831 //=======================================================================
2832 ///Shape functions for PRefineableQElement<DIM>
2833 //=======================================================================
2834 template<unsigned INITIAL_NNODE_1D>
2836 shape(const Vector<double> &s, Shape &psi) const
2837 {
2838  switch(p_order())
2839  {
2840  case 2:
2841  {
2842  //Call the OneDimensional Shape functions
2844  OneDimensionalLegendreShape<2> psi1(s[0]), psi2(s[1]);
2845 
2846  //Now let's loop over the nodal points in the element
2847  //and copy the values back in
2848  for(unsigned i=0;i<2;i++)
2849  {
2850  for(unsigned j=0;j<2;j++)
2851  {
2852  psi(2*i + j) = psi2[i]*psi1[j];
2853  }
2854  }
2855  break;
2856  }
2857  case 3:
2858  {
2859  //Call the OneDimensional Shape functions
2861  OneDimensionalLegendreShape<3> psi1(s[0]), psi2(s[1]);
2862 
2863  //Now let's loop over the nodal points in the element
2864  //and copy the values back in
2865  for(unsigned i=0;i<3;i++)
2866  {
2867  for(unsigned j=0;j<3;j++)
2868  {
2869  psi(3*i + j) = psi2[i]*psi1[j];
2870  }
2871  }
2872  break;
2873  }
2874  case 4:
2875  {
2876  //Call the OneDimensional Shape functions
2878  OneDimensionalLegendreShape<4> psi1(s[0]), psi2(s[1]);
2879 
2880  //Now let's loop over the nodal points in the element
2881  //and copy the values back in
2882  for(unsigned i=0;i<4;i++)
2883  {
2884  for(unsigned j=0;j<4;j++)
2885  {
2886  psi(4*i + j) = psi2[i]*psi1[j];
2887  }
2888  }
2889  break;
2890  }
2891  case 5:
2892  {
2893  //Call the OneDimensional Shape functions
2895  OneDimensionalLegendreShape<5> psi1(s[0]), psi2(s[1]);
2896 
2897  //Now let's loop over the nodal points in the element
2898  //and copy the values back in
2899  for(unsigned i=0;i<5;i++)
2900  {
2901  for(unsigned j=0;j<5;j++)
2902  {
2903  psi(5*i + j) = psi2[i]*psi1[j];
2904  }
2905  }
2906  break;
2907  }
2908  case 6:
2909  {
2910  //Call the OneDimensional Shape functions
2912  OneDimensionalLegendreShape<6> psi1(s[0]), psi2(s[1]);
2913 
2914  //Now let's loop over the nodal points in the element
2915  //and copy the values back in
2916  for(unsigned i=0;i<6;i++)
2917  {
2918  for(unsigned j=0;j<6;j++)
2919  {
2920  psi(6*i + j) = psi2[i]*psi1[j];
2921  }
2922  }
2923  break;
2924  }
2925  case 7:
2926  {
2927  //Call the OneDimensional Shape functions
2929  OneDimensionalLegendreShape<7> psi1(s[0]), psi2(s[1]);
2930 
2931  //Now let's loop over the nodal points in the element
2932  //and copy the values back in
2933  for(unsigned i=0;i<7;i++)
2934  {
2935  for(unsigned j=0;j<7;j++)
2936  {
2937  psi(7*i + j) = psi2[i]*psi1[j];
2938  }
2939  }
2940  break;
2941  }
2942  default:
2943  std::ostringstream error_message;
2944  error_message <<"\nERROR: Exceeded maximum polynomial order for";
2945  error_message <<"\n polynomial order for shape functions.\n";
2946  throw OomphLibError(error_message.str(),
2947  OOMPH_CURRENT_FUNCTION,
2948  OOMPH_EXCEPTION_LOCATION);
2949  }
2950 
2951 }
2952 
2953 //=======================================================================
2954 ///Derivatives of shape functions for PRefineableQElement<DIM>
2955 //=======================================================================
2956 template<unsigned INITIAL_NNODE_1D>
2958 dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsids) const
2959 {
2960  switch(p_order())
2961  {
2962  case 2:
2963  {
2964  //Call the shape functions and derivatives
2966  OneDimensionalLegendreShape<2> psi1(s[0]), psi2(s[1]);
2967  OneDimensionalLegendreDShape<2> dpsi1ds(s[0]), dpsi2ds(s[1]);
2968 
2969  //Index for the shape functions
2970  unsigned index=0;
2971  //Loop over shape functions in element
2972  for(unsigned i=0;i<2;i++)
2973  {
2974  for(unsigned j=0;j<2;j++)
2975  {
2976  //Assign the values
2977  dpsids(index,0) = psi2[i]*dpsi1ds[j];
2978  dpsids(index,1) = dpsi2ds[i]*psi1[j];
2979  psi[index] = psi2[i]*psi1[j];
2980  //Increment the index
2981  ++index;
2982  }
2983  }
2984  break;
2985  }
2986  case 3:
2987  {
2988  //Call the shape functions and derivatives
2990  OneDimensionalLegendreShape<3> psi1(s[0]), psi2(s[1]);
2991  OneDimensionalLegendreDShape<3> dpsi1ds(s[0]), dpsi2ds(s[1]);
2992 
2993  //Index for the shape functions
2994  unsigned index=0;
2995  //Loop over shape functions in element
2996  for(unsigned i=0;i<3;i++)
2997  {
2998  for(unsigned j=0;j<3;j++)
2999  {
3000  //Assign the values
3001  dpsids(index,0) = psi2[i]*dpsi1ds[j];
3002  dpsids(index,1) = dpsi2ds[i]*psi1[j];
3003  psi[index] = psi2[i]*psi1[j];
3004  //Increment the index
3005  ++index;
3006  }
3007  }
3008  break;
3009  }
3010  case 4:
3011  {
3012  //Call the shape functions and derivatives
3014  OneDimensionalLegendreShape<4> psi1(s[0]), psi2(s[1]);
3015  OneDimensionalLegendreDShape<4> dpsi1ds(s[0]), dpsi2ds(s[1]);
3016 
3017  //Index for the shape functions
3018  unsigned index=0;
3019  //Loop over shape functions in element
3020  for(unsigned i=0;i<4;i++)
3021  {
3022  for(unsigned j=0;j<4;j++)
3023  {
3024  //Assign the values
3025  dpsids(index,0) = psi2[i]*dpsi1ds[j];
3026  dpsids(index,1) = dpsi2ds[i]*psi1[j];
3027  psi[index] = psi2[i]*psi1[j];
3028  //Increment the index
3029  ++index;
3030  }
3031  }
3032  break;
3033  }
3034  case 5:
3035  {
3036  //Call the shape functions and derivatives
3038  OneDimensionalLegendreShape<5> psi1(s[0]), psi2(s[1]);
3039  OneDimensionalLegendreDShape<5> dpsi1ds(s[0]), dpsi2ds(s[1]);
3040 
3041  //Index for the shape functions
3042  unsigned index=0;
3043  //Loop over shape functions in element
3044  for(unsigned i=0;i<5;i++)
3045  {
3046  for(unsigned j=0;j<5;j++)
3047  {
3048  //Assign the values
3049  dpsids(index,0) = psi2[i]*dpsi1ds[j];
3050  dpsids(index,1) = dpsi2ds[i]*psi1[j];
3051  psi[index] = psi2[i]*psi1[j];
3052  //Increment the index
3053  ++index;
3054  }
3055  }
3056  break;
3057  }
3058  case 6:
3059  {
3060  //Call the shape functions and derivatives
3062  OneDimensionalLegendreShape<6> psi1(s[0]), psi2(s[1]);
3063  OneDimensionalLegendreDShape<6> dpsi1ds(s[0]), dpsi2ds(s[1]);
3064 
3065  //Index for the shape functions
3066  unsigned index=0;
3067  //Loop over shape functions in element
3068  for(unsigned i=0;i<6;i++)
3069  {
3070  for(unsigned j=0;j<6;j++)
3071  {
3072  //Assign the values
3073  dpsids(index,0) = psi2[i]*dpsi1ds[j];
3074  dpsids(index,1) = dpsi2ds[i]*psi1[j];
3075  psi[index] = psi2[i]*psi1[j];
3076  //Increment the index
3077  ++index;
3078  }
3079  }
3080  break;
3081  }
3082  case 7:
3083  {
3084  //Call the shape functions and derivatives
3086  OneDimensionalLegendreShape<7> psi1(s[0]), psi2(s[1]);
3087  OneDimensionalLegendreDShape<7> dpsi1ds(s[0]), dpsi2ds(s[1]);
3088 
3089  //Index for the shape functions
3090  unsigned index=0;
3091  //Loop over shape functions in element
3092  for(unsigned i=0;i<7;i++)
3093  {
3094  for(unsigned j=0;j<7;j++)
3095  {
3096  //Assign the values
3097  dpsids(index,0) = psi2[i]*dpsi1ds[j];
3098  dpsids(index,1) = dpsi2ds[i]*psi1[j];
3099  psi[index] = psi2[i]*psi1[j];
3100  //Increment the index
3101  ++index;
3102  }
3103  }
3104  break;
3105  }
3106  default:
3107  std::ostringstream error_message;
3108  error_message <<"\nERROR: Exceeded maximum polynomial order for";
3109  error_message <<"\n polynomial order for shape functions.\n";
3110  throw OomphLibError(error_message.str(),
3111  OOMPH_CURRENT_FUNCTION,
3112  OOMPH_EXCEPTION_LOCATION);
3113  }
3114 
3115 }
3116 
3117 //=======================================================================
3118 /// Second derivatives of shape functions for PRefineableQElement<DIM>
3119 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
3120 //=======================================================================
3121 template<unsigned INITIAL_NNODE_1D>
3124  DShape &d2psids) const
3125 {
3126  std::ostringstream error_message;
3127  error_message <<"\nd2shape_local currently not implemented for this element\n";
3128  throw OomphLibError(error_message.str(),
3129  OOMPH_CURRENT_FUNCTION,
3130  OOMPH_EXCEPTION_LOCATION);
3131 }
3132 
3133 //=======================================================================
3134 /// Rebuild the element from nodes found in its sons
3135 /// Adjusts its p-order to be the maximum of its sons' p-orders
3136 //=======================================================================
3137 template<unsigned INITIAL_NNODE_1D>
3140 {
3141  using namespace QuadTreeNames;
3142 
3143  // Get p-orders of sons
3144  unsigned n_sons = this->tree_pt()->nsons();
3145  Vector<unsigned> son_p_order(n_sons);
3146  unsigned max_son_p_order = 0;
3147  for (unsigned ison=0;ison<n_sons;ison++)
3148  {
3149  PRefineableElement* el_pt = dynamic_cast<PRefineableElement*>(this->tree_pt()->son_pt(ison)->object_pt());
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];
3152  }
3153 
3154  unsigned old_Nnode = this->nnode();
3155  unsigned old_P_order = this->p_order();
3156  // Set p-order of the element
3157  this->p_order() = max_son_p_order;
3158 
3159  // Change integration scheme
3160  delete this->integral_pt();
3161  switch(this->p_order())
3162  {
3163  case 2:
3164  this->set_integration_scheme(new GaussLobattoLegendre<2,2>);
3165  break;
3166  case 3:
3167  this->set_integration_scheme(new GaussLobattoLegendre<2,3>);
3168  break;
3169  case 4:
3170  this->set_integration_scheme(new GaussLobattoLegendre<2,4>);
3171  break;
3172  case 5:
3173  this->set_integration_scheme(new GaussLobattoLegendre<2,5>);
3174  break;
3175  case 6:
3176  this->set_integration_scheme(new GaussLobattoLegendre<2,6>);
3177  break;
3178  case 7:
3179  this->set_integration_scheme(new GaussLobattoLegendre<2,7>);
3180  break;
3181  default:
3182  std::ostringstream error_message;
3183  error_message <<"\nERROR: Exceeded maximum polynomial order for";
3184  error_message <<"\n integration scheme.\n";
3185  throw OomphLibError(error_message.str(),
3186  OOMPH_CURRENT_FUNCTION,
3187  OOMPH_EXCEPTION_LOCATION);
3188  }
3189 
3190  // Back up pointers to old nodes before they are lost
3191  Vector<Node*> old_node_pt(old_Nnode);
3192  for (unsigned n=0; n<old_Nnode; n++)
3193  {
3194  old_node_pt[n] = this->node_pt(n);
3195  }
3196 
3197  // Allocate new space for Nodes (at the element level)
3198  this->set_n_node(this->p_order()*this->p_order());
3199 
3200  // Copy vertex nodes which were populated in the pre-build
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];
3207 
3208  // Copy midpoint nodes from sons if new p-order is odd
3209  if(this->p_order() % 2 == 1)
3210  {
3211  //Work out which is midpoint node
3212  unsigned midpoint = (this->p_order()-1)/2;
3213 
3214  //Bottom edge
3215  this->node_pt(midpoint)
3216  = dynamic_cast<RefineableQElement<2>*>
3217  (quadtree_pt()->son_pt(SW)->object_pt())->vertex_node_pt(1);
3218  //Left edge
3219  this->node_pt(midpoint*this->p_order())
3220  = dynamic_cast<RefineableQElement<2>*>
3221  (quadtree_pt()->son_pt(SW)->object_pt())->vertex_node_pt(2);
3222  //Top edge
3223  this->node_pt((this->p_order()-1)*this->p_order()+midpoint)
3224  = dynamic_cast<RefineableQElement<2>*>
3225  (quadtree_pt()->son_pt(NE)->object_pt())->vertex_node_pt(2);
3226  //Right edge
3227  this->node_pt((midpoint+1)*this->p_order()-1)
3228  = dynamic_cast<RefineableQElement<2>*>
3229  (quadtree_pt()->son_pt(NE)->object_pt())->vertex_node_pt(1);
3230  }
3231 
3232 
3233 
3234 
3235  //The timestepper should be the same for all nodes and node 0 should
3236  //never be deleted.
3237  if(this->node_pt(0)==0)
3238  {
3239  throw OomphLibError("The Corner node (0) does not exist",
3240  OOMPH_CURRENT_FUNCTION,
3241  OOMPH_EXCEPTION_LOCATION);
3242  }
3243 
3244  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
3245  unsigned ntstorage = time_stepper_pt->ntstorage();
3246 
3247  unsigned jnod=0;
3248  Vector<double> s_fraction(2), s(2);
3249  //Loop over the nodes in the element
3250  unsigned n_p = this->nnode_1d();
3251  for(unsigned i0=0;i0<n_p;i0++)
3252  {
3253  //Get the fractional position of the node
3254  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
3255  //Local coordinate
3256  s[0] = -1.0 + 2.0*s_fraction[0];
3257 
3258  for(unsigned i1=0;i1<n_p;i1++)
3259  {
3260  //Get the fractional position of the node in the direction of s[1]
3261  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
3262  // Local coordinate in father element
3263  s[1] = -1.0 + 2.0*s_fraction[1];
3264 
3265  //Set the local node number
3266  jnod = i0 + n_p*i1;
3267 
3268  // Initialise flag: So far, this node hasn't been built
3269  // or copied yet
3270  bool node_done=false;
3271 
3272  //Get the pointer to the node in this element, returns NULL
3273  //if there is not node
3274  Node* created_node_pt = this->get_node_at_local_coordinate(s);
3275 
3276  // Does this node already exist in this element?
3277  //----------------------------------------------
3278  if (created_node_pt!=0)
3279  {
3280  // Copy node across
3281  this->node_pt(jnod) = created_node_pt;
3282 
3283  // Node has been created by copy
3284  node_done = true;
3285  }
3286  // Node does not exist in this element but might already
3287  //------------------------------------------------------
3288  // have been created by neighbouring elements
3289  //-------------------------------------------
3290  else
3291  {
3292  //Was the node created by one of its neighbours
3293  //Whether or not the node lies on an edge can be calculated
3294  //by from the fractional position
3295  bool is_periodic = false;
3296  created_node_pt =
3297  node_created_by_neighbour(s_fraction,is_periodic);
3298 
3299  //If the node was so created, assign the pointers
3300  if(created_node_pt!=0)
3301  {
3302  //If the node is periodic
3303  if(is_periodic)
3304  {
3305  throw OomphLibError(
3306  "Cannot handle periodic nodes yet",
3307  OOMPH_CURRENT_FUNCTION,
3308  OOMPH_EXCEPTION_LOCATION);
3309  }
3310  //Non-periodic case, just set the pointer
3311  else
3312  {
3313  this->node_pt(jnod) = created_node_pt;
3314  }
3315  //Node has been created
3316  node_done = true;
3317  }
3318  } // Node does not exist in this element
3319 
3320  // Node has not been built anywhere ---> build it here
3321  if (!node_done)
3322  {
3323  //First, find the son element in which the node should live
3324 
3325  //Find coordinates in the sons
3326  Vector<double> s_in_son(2);
3327  using namespace QuadTreeNames;
3328  int son=-10;
3329  //If negative on the west side
3330  if(s_fraction[0] < 0.5)
3331  {
3332  //On the south side
3333  if(s_fraction[1] < 0.5)
3334  {
3335  //It's the southwest son
3336  son = SW;
3337  s_in_son[0] = -1.0 + 4.0*s_fraction[0];
3338  s_in_son[1] = -1.0 + 4.0*s_fraction[1];
3339  }
3340  //On the north side
3341  else
3342  {
3343  //It's the northwest son
3344  son = NW;
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);
3347  }
3348  }
3349  else
3350  {
3351  //On the south side
3352  if(s_fraction[1] < 0.5)
3353  {
3354  //It's the southeast son
3355  son = SE;
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];
3358  }
3359  //On the north side
3360  else
3361  {
3362  //It's the northeast son
3363  son = NE;
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);
3366  }
3367  }
3368 
3369  //Get the pointer to the son element in which the new node
3370  //would live
3373  this->tree_pt()->son_pt(son)->object_pt());
3374 
3375  //If we are rebuilding, then worry about the boundary conditions
3376  //Find the boundary of the node
3377  //Initially none
3378  int boundary=Tree::OMEGA;
3379  //If we are on the western boundary
3380  if(i0==0) {boundary = W;}
3381  //If we are on the eastern boundary
3382  else if(i0==n_p-1) {boundary = E;}
3383 
3384  //If we are on the southern boundary
3385  if(i1==0)
3386  {
3387  //If we already have already set the boundary, we're on a corner
3388  switch(boundary)
3389  {
3390  case W:
3391  boundary = SW;
3392  break;
3393  case E:
3394  boundary = SE;
3395  break;
3396  //Boundary not set
3397  default:
3398  boundary = S;
3399  break;
3400  }
3401  }
3402  //If we are the northern bounadry
3403  else if(i1==n_p-1)
3404  {
3405  //If we already have a boundary
3406  switch(boundary)
3407  {
3408  case W:
3409  boundary = NW;
3410  break;
3411  case E:
3412  boundary = NE;
3413  break;
3414  default:
3415  boundary = N;
3416  break;
3417  }
3418  }
3419 
3420  // set of boundaries that this edge in the son lives on
3421  std::set<unsigned> boundaries;
3422 
3423  //Now get the boundary conditions from the son
3424  //The boundaries will be common to the son because there can be
3425  //no rotations here
3426  if(boundary!=Tree::OMEGA)
3427  {
3428  son_el_pt->get_boundaries(boundary,boundaries);
3429  }
3430 
3431  // If the node lives on a boundary:
3432  // Construct a boundary node,
3433  // Get boundary conditions and
3434  // update all lookup schemes
3435  if(boundaries.size()>0)
3436  {
3437  //Construct the new node
3438  created_node_pt = this->construct_boundary_node(jnod,time_stepper_pt);
3439 
3440  //Get the boundary conditions from the son
3441  Vector<int> bound_cons(this->ncont_interpolated_values());
3442  son_el_pt->get_bcs(boundary,bound_cons);
3443 
3444  //Loop over the values and pin if necessary
3445  unsigned nval = created_node_pt->nvalue();
3446  for(unsigned k=0;k<nval;k++)
3447  {
3448  if(bound_cons[k]) {created_node_pt->pin(k);}
3449  }
3450 
3451  // Solid node? If so, deal with the positional boundary
3452  // conditions:
3453  SolidNode* solid_node_pt =
3454  dynamic_cast<SolidNode*>(created_node_pt);
3455  if (solid_node_pt!=0)
3456  {
3457  //Get the positional boundary conditions from the father:
3458  unsigned n_dim = created_node_pt->ndim();
3459  Vector<int> solid_bound_cons(n_dim);
3460  RefineableSolidQElement<2>* son_solid_el_pt=
3461  dynamic_cast<RefineableSolidQElement<2>*>(son_el_pt);
3462 #ifdef PARANOID
3463  if (son_solid_el_pt==0)
3464  {
3465  std::string error_message =
3466  "We have a SolidNode outside a refineable SolidElement\n";
3467  error_message +=
3468  "during mesh refinement -- this doesn't make sense\n";
3469 
3470  throw OomphLibError(
3471  error_message,
3472  OOMPH_CURRENT_FUNCTION,
3473  OOMPH_EXCEPTION_LOCATION);
3474  }
3475 #endif
3476  son_solid_el_pt->
3477  get_solid_bcs(boundary,solid_bound_cons);
3478 
3479  //Loop over the positions and pin, if necessary
3480  for(unsigned k=0;k<n_dim;k++)
3481  {
3482  if (solid_bound_cons[k]) {solid_node_pt->pin_position(k);}
3483  }
3484  } //End of if solid_node_pt
3485 
3486 
3487 
3488  //Next we update the boundary look-up schemes
3489  //Loop over the boundaries stored in the set
3490  for(std::set<unsigned>::iterator it = boundaries.begin();
3491  it != boundaries.end(); ++it)
3492  {
3493  //Add the node to the boundary
3494  mesh_pt->add_boundary_node(*it,created_node_pt);
3495 
3496  //If we have set an intrinsic coordinate on this
3497  //mesh boundary then it must also be interpolated on
3498  //the new node
3499  //Now interpolate the intrinsic boundary coordinate
3500  if(mesh_pt->boundary_coordinate_exists(*it)==true)
3501  {
3502  Vector<double> zeta(1);
3503  son_el_pt->interpolated_zeta_on_edge(*it,boundary,
3504  s_in_son,zeta);
3505 
3506  created_node_pt->set_coordinates_on_boundary(*it,zeta);
3507  }
3508  }
3509  }
3510  //Otherwise the node is not on a Mesh boundary
3511  //and we create a normal "bulk" node
3512  else
3513  {
3514  //Construct the new node
3515  created_node_pt = this->construct_node(jnod,time_stepper_pt);
3516  }
3517 
3518  //Now we set the position and values at the newly created node
3519 
3520  // In the first instance use macro element or FE representation
3521  // to create past and present nodal positions.
3522  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
3523  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
3524  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
3525  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
3526  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
3527  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
3528  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
3529 
3530  // Loop over # of history values
3531  //Loop over # of history values
3532  for(unsigned t=0;t<ntstorage;t++)
3533  {
3534  using namespace QuadTreeNames;
3535  //Get the position from the son
3536  Vector<double> x_prev(2);
3537 
3538  //Now let's fill in the value
3539  son_el_pt->get_x(t,s_in_son,x_prev);
3540  for(unsigned i=0;i<2;i++)
3541  {
3542  created_node_pt->x(t,i) = x_prev[i];
3543  }
3544  }
3545 
3546  // Now set up the values
3547  // Loop over all history values
3548  for(unsigned t=0;t<ntstorage;t++)
3549  {
3550  // Get values from father element
3551  // Note: get_interpolated_values() sets Vector size itself.
3552  Vector<double> prev_values;
3553  son_el_pt->get_interpolated_values(t,s_in_son,prev_values);
3554 
3555  //Initialise the values at the new node
3556  for(unsigned k=0;k<created_node_pt->nvalue();k++)
3557  {
3558  created_node_pt->set_value(t,k,prev_values[k]);
3559  }
3560  }
3561 
3562  //Add the node to the mesh
3563  mesh_pt->add_node_pt(created_node_pt);
3564 
3565  // Check if the element is an algebraic element
3566  AlgebraicElementBase* alg_el_pt =
3567  dynamic_cast<AlgebraicElementBase*>(this);
3568 
3569  //If we do have an algebraic element
3570  if(alg_el_pt!=0)
3571  {
3572  std::string error_message =
3573  "Have not implemented rebuilding from sons for";
3574  error_message +=
3575  "Algebraic p-refineable elements yet\n";
3576 
3577  throw
3578  OomphLibError(error_message,
3579  "PRefineableQElement<2,INITIAL_NNODE_1D>::rebuild_from_sons()",
3580  OOMPH_EXCEPTION_LOCATION);
3581  }
3582 
3583  } //End of the case when we build the node ourselves
3584  }
3585  } // End of loop over all nodes in element
3586 
3587 
3588  // If the element is a MacroElementNodeUpdateElement, set the update
3589  // parameters for the current element's nodes. These need to be reset
3590  // (as in MacroElementNodeUpdateElement<ELEMENT>::rebuild_from_sons())
3591  // because the nodes in this element have changed
3592  MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
3594  if(m_el_pt!=0)
3595  {
3596  // Loop over the nodes
3597  for (unsigned j=0;j<this->nnode();j++)
3598  {
3599  // Get local coordinate in element (Vector sets its own size)
3600  Vector<double> s_in_node_update_element;
3601  this->local_coordinate_of_node(j,s_in_node_update_element);
3602 
3603  // Get vector of geometric objects
3604  Vector<GeomObject*> geom_object_pt(m_el_pt->geom_object_pt());
3605 
3606  // Pass the lot to the node
3607  static_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))->
3608  set_node_update_info(this,s_in_node_update_element,geom_object_pt);
3609  }
3610  }
3611 
3612  //MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
3613  // MacroElementNodeUpdateElementBase*>(this);
3614  //if(m_el_pt!=0)
3615  // {
3616  // Loop over all nodes in element again, to re-set the positions
3617  // This must be done using the new element's macro-element
3618  // representation, rather than the old version which may be
3619  // of a different p-order!
3620  for(unsigned i0=0;i0<n_p;i0++)
3621  {
3622  //Get the fractional position of the node
3623  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
3624  //Local coordinate
3625  s[0] = -1.0 + 2.0*s_fraction[0];
3626 
3627  for(unsigned i1=0;i1<n_p;i1++)
3628  {
3629  //Get the fractional position of the node in the direction of s[1]
3630  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
3631  // Local coordinate in father element
3632  s[1] = -1.0 + 2.0*s_fraction[1];
3633 
3634  //Set the local node number
3635  jnod = i0 + n_p*i1;
3636 
3637  // Loop over # of history values
3638  for (unsigned t=0;t<ntstorage;t++)
3639  {
3640  // Get position from father element -- this uses the macro
3641  // element representation if appropriate. If the node
3642  // turns out to be a hanging node later on, then
3643  // its position gets adjusted in line with its
3644  // hanging node interpolation.
3645  Vector<double> x_prev(2);
3646  this->get_x(t,s,x_prev);
3647 
3648  // Set previous positions of the new node
3649  for(unsigned i=0;i<2;i++)
3650  {
3651  this->node_pt(jnod)->x(t,i) = x_prev[i];
3652  }
3653  }
3654  }
3655  }
3656  // }
3657 
3658 }
3659 
3660 //=================================================================
3661 /// Check inter-element continuity of
3662 /// - nodal positions
3663 /// - (nodally) interpolated function values
3664 /// Overloaded to not check differences in the value. Mortaring
3665 /// doesn't enforce strong continuity between elements.
3666 //====================================================================
3667 template<unsigned INITIAL_NNODE_1D>
3669 check_integrity(double& max_error)
3670 {
3671  // Overloaded to *not* check for continuity in value of interpolated
3672  // variables. This is necessary because mortaring does not ensure continuity
3673  // across element boundaries. It therefore makes no sense to test for this.
3674 
3675  //Dummy set max_error to 0
3676  max_error = 0.0;
3677 
3678  // With macro-elements, (strong) continuity in position is nolonger
3679  // guaranteed either, so we don't check for this either. In fact, we do
3680  // nothing at all.
3681  if(this->macro_elem_pt()!=0)
3682  {
3683  //We have a macro element, so do nothing!
3684  return;
3685  }
3686 
3687  using namespace QuadTreeNames;
3688 
3689  // Number of nodes along edge
3690  unsigned n_p=nnode_1d();
3691 
3692  // Number of timesteps (incl. present) for which continuity is
3693  // to be checked.
3694  unsigned n_time=1;
3695 
3696  // Initialise errors
3697  max_error=0.0;
3698  Vector<double> max_error_x(2,0.0);
3699  double max_error_val=0.0;
3700 
3701  Vector<int> edges(4);
3702  edges[0] = S; edges[1] = N; edges[2] = W; edges[3] = E;
3703 
3704  //Loop over the edges
3705  for(unsigned edge_counter=0;edge_counter<4;edge_counter++)
3706  {
3707  Vector<unsigned> translate_s(2);
3708  Vector<double> s(2), s_lo_neigh(2), s_hi_neigh(2), s_fraction(2);
3709  int neigh_edge,diff_level;
3710  bool in_neighbouring_tree;
3711 
3712  // Find pointer to neighbour in this direction
3713  QuadTree* neigh_pt;
3714  neigh_pt=quadtree_pt()->gteq_edge_neighbour(edges[edge_counter],
3715  translate_s,
3716  s_lo_neigh,s_hi_neigh,
3717  neigh_edge,diff_level,
3718  in_neighbouring_tree);
3719 
3720  // Neighbour exists and has existing nodes
3721  if((neigh_pt!=0) && (neigh_pt->object_pt()->nodes_built()))
3722  {
3723  //Need to exclude periodic nodes from this check
3724  //There are only periodic nodes if we are in a neighbouring tree
3725  bool is_periodic=false;
3726  if(in_neighbouring_tree)
3727  {
3728  //Is it periodic
3729  is_periodic =
3730  this->tree_pt()->root_pt()->is_neighbour_periodic(edges[edge_counter]);
3731  }
3732 
3733  // We also need to exclude edges which may have hanging nodes
3734  // because mortaring does not guarantee (strong) continuity
3735  // in position or in value at nonconforming element boundaries
3736  bool exclude_this_edge = false;
3737  if(diff_level != 0)
3738  {
3739  // h-type nonconformity (slave)
3740  exclude_this_edge = true;
3741  }
3742  else if(neigh_pt->nsons() != 0)
3743  {
3744  // h-type nonconformity (master)
3745  exclude_this_edge = true;
3746  }
3747  else
3748  {
3749  unsigned my_p_order = this->p_order();
3750  unsigned neigh_p_order =
3751  dynamic_cast<PRefineableQElement*>(neigh_pt->object_pt())->p_order();
3752  if(my_p_order != neigh_p_order)
3753  {
3754  // p-type nonconformity
3755  exclude_this_edge = true;
3756  }
3757  }
3758 
3759  // With macro-elements, (strong) continuity in position is nolonger
3760  // guaranteed either, so we don't check for this either. In fact, we do
3761  // nothing at all.
3762  if(dynamic_cast<FiniteElement*>
3763  (neigh_pt->object_pt())->macro_elem_pt()!=0)
3764  {
3765  //We have a macro element, so do nothing!
3766  break;
3767  }
3768 
3769  //Check conforming edges
3770  if(!exclude_this_edge)
3771  {
3772 
3773  // Loop over nodes along the edge
3774  for(unsigned i0=0;i0<n_p;i0++)
3775  {
3776  //Storage for pointer to the local node
3777  Node* local_node_pt=0;
3778 
3779  switch(edge_counter)
3780  {
3781  case 0:
3782  // Local fraction of node
3783  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
3784  s_fraction[1] = 0.0;
3785  // Get pointer to local node
3786  local_node_pt = this->node_pt(i0);
3787  break;
3788 
3789  case 1:
3790  // Local fraction of node
3791  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
3792  s_fraction[1] = 1.0;
3793  // Get pointer to local node
3794  local_node_pt = this->node_pt(i0 + n_p*(n_p-1));
3795  break;
3796 
3797  case 2:
3798  // Local fraction of node
3799  s_fraction[0] = 0.0;
3800  s_fraction[1] = this->local_one_d_fraction_of_node(i0,1);
3801  // Get pointer to local node
3802  local_node_pt = this->node_pt(n_p*i0);
3803  break;
3804 
3805  case 3:
3806  // Local fraction of node
3807  s_fraction[0] = 1.0;
3808  s_fraction[1] = this->local_one_d_fraction_of_node(i0,1);
3809  // Get pointer to local node
3810  local_node_pt = this->node_pt(n_p-1 + n_p*i0);
3811  break;
3812  }
3813 
3814  //Calculate the local coordinate and the local coordinate as viewed
3815  //from the neighbour
3816  Vector<double> s_in_neighb(2);
3817  for(unsigned i=0;i<2;i++)
3818  {
3819  //Local coordinate in this element
3820  s[i] = -1.0 + 2.0*s_fraction[i];
3821  //Local coordinate in the neighbour
3822  s_in_neighb[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
3823  (s_hi_neigh[i] - s_lo_neigh[i]);
3824  }
3825 
3826  //Loop over timesteps
3827  for(unsigned t=0;t<n_time;t++)
3828  {
3829  // Get the nodal position from neighbour element
3830  Vector<double> x_in_neighb(2);
3831  neigh_pt->object_pt()->interpolated_x(t,s_in_neighb,x_in_neighb);
3832 
3833  // Check error only if the node is NOT periodic
3834  if(is_periodic==false)
3835  {
3836  for(int i=0;i<2;i++)
3837  {
3838  //Find the spatial error
3839  double err = std::fabs(local_node_pt->x(t,i) - x_in_neighb[i]);
3840 
3841  //If it's bigger than our tolerance, say so
3842  if (err>1e-9)
3843  {
3844  oomph_info << "errx " << err << " " << t << " "
3845  << local_node_pt->x(t,i)
3846  << " " << x_in_neighb[i]<< std::endl;
3847 
3848  oomph_info << "at " << local_node_pt->x(0) << " "
3849  << local_node_pt->x(1) << std::endl;
3850  }
3851 
3852  //If it's bigger than the previous max error, it is the
3853  //new max error!
3854  if (err>max_error_x[i]) {max_error_x[i]=err;}
3855  }
3856  }
3857 
3858  // Get the values from neighbour element. Note: # of values
3859  // gets set by routine (because in general we don't know
3860  // how many interpolated values a certain element has
3861  Vector<double> values_in_neighb;
3862  neigh_pt->object_pt()->
3863  get_interpolated_values(t,s_in_neighb,values_in_neighb);
3864 
3865  // Get the values in current element.
3866  Vector<double> values;
3867  this->get_interpolated_values(t,s,values);
3868 
3869  // Now figure out how many continuously interpolated values there are
3870  unsigned num_val=neigh_pt->object_pt()->ncont_interpolated_values();
3871 
3872  // Check error
3873  for(unsigned ival=0;ival<num_val;ival++)
3874  {
3875  double err=std::fabs(values[ival] - values_in_neighb[ival]);
3876 
3877  if (err>1.0e-10)
3878  {
3879  oomph_info << local_node_pt->x(0) << " "
3880  << local_node_pt->x(1) << " \n# "
3881  << "erru (S)" << err << " " << ival << " "
3882  << this->get_node_number(local_node_pt) << " "
3883  << values[ival]
3884  << " " << values_in_neighb[ival] << std::endl;
3885  }
3886 
3887  if (err>max_error_val) {max_error_val=err;}
3888 
3889  }
3890  }
3891 
3892  }
3893  }
3894  }
3895  }
3896 
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;
3900 
3901  if (max_error>1e-9)
3902  {
3903  oomph_info << "\n#------------------------------------ \n#Max error " ;
3904  oomph_info << max_error_x[0]
3905  << " " << max_error_x[1]
3906  << " " << max_error_val << std::endl;
3907  oomph_info << "#------------------------------------ \n " << std::endl;
3908 
3909  }
3910 
3911 }
3912 
3913 //=================================================================
3914 /// Internal function to set up the hanging nodes on a particular
3915 /// edge of the element.
3916 /// Implements the mortarting method to enforce continuity weakly
3917 /// across non-conforming element boundaries \f$\Gamma\f$ using an
3918 /// integral matching condition
3919 /// \f[ \int_\Gamma (u_{\mbox{S}} - u_{\mbox{M}}) \psi \mbox{d} s = 0 \f]
3920 /// for all polynomials \f$\psi\f$ on \f$\Gamma\f$ of degree at most
3921 /// p-2 (where p is the spectral-order of the slave element) and a
3922 /// vertex matching condition
3923 /// \f[ (u_{\mbox{S}} - u_{\mbox{M}})\big\vert_{\partial\Gamma} = 0.\f]
3924 ///
3925 /// The algorithm works as follows:
3926 /// - First the element determines if its edge my_edge is on the
3927 /// master or slave side of the non-conformity. At h-type non-conformities
3928 /// we choose long edges to be masters, and at p-type nonconformities the
3929 /// edge with lower p-order is the master.
3930 /// - Mortaring is performed by the slave side.
3931 /// - If a vertex node of the mortar is shared between master and slave
3932 /// element then the vertex matching condition is enforced automatically.
3933 /// Otherwise it must be imposed by constraining its value to that at on
3934 /// the master side.
3935 /// - The integral matching condition is discretised and the mortar test
3936 /// functions \f$ \psi \f$ are chosen to be derivatives of Legendre
3937 /// polynomials of degree p-1.
3938 /// - The mortar mass matrix M is constructed. Its entries are the
3939 /// mortar test functions evaluated at the slave nodal positions, so it is
3940 /// diagonal.
3941 /// - The local projection matrix is constructed for the master element by
3942 /// applying the discretised integral matching condition along the mortar
3943 /// using the appropriate quadrature order.
3944 /// - The global projection matrix is then assembled by subtracting
3945 /// contributions from the mortar vertex nodes.
3946 /// - The mortar system \f$ M\xi^s = P\hat{\xi^m} \f$ is constructed,
3947 /// where \f$ \xi^m \f$ and \f$ \xi^s \f$ are the nodal values at the master
3948 /// and slave nodes respectively.
3949 /// - The conformity matrix \f$ C = M^{-1}P \f$ is computed. This is
3950 /// straightforward since the mass matrix is diagonal.
3951 /// - Finally, the master nodes and weights for each slave node are read from
3952 /// the conformity matrix and stored in the slaves' hanging schemes.
3953 ///
3954 /// The positions of the slave nodes are set to be consistent with their
3955 /// hanging schemes.
3956 //=================================================================
3957 template<unsigned INITIAL_NNODE_1D>
3959 quad_hang_helper(const int &value_id, const int &my_edge,
3960  std::ofstream& output_hangfile)
3961 {
3962  using namespace QuadTreeNames;
3963 
3964  Vector<unsigned> translate_s(2);
3965  Vector<double> s_lo_neigh(2);
3966  Vector<double> s_hi_neigh(2);
3967  int neigh_edge,diff_level;
3968  bool in_neighbouring_tree;
3969 
3970  // Find pointer to neighbour in this direction
3971  QuadTree* neigh_pt;
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);
3975 
3976  // Work out master/slave edges
3977  //----------------------------
3978 
3979  // Set up booleans
3980  //bool h_type_master = false;
3981  bool h_type_slave = false;
3982  //bool p_type_master = false;
3983  bool p_type_slave = false;
3984 
3985  // Neighbour exists and all nodes have been created
3986  if(neigh_pt!=0)
3987  {
3988  // Check if neighbour is bigger than me
3989  if(diff_level!=0)
3990  {
3991  // Slave at h-type non-conformity
3992  h_type_slave = true;
3993  }
3994  // Check if neighbour is the same size as me
3995  else if(neigh_pt->nsons()==0)
3996  {
3997  // Neighbour is same size as me
3998  // Find p-orders of each element
3999  unsigned my_p_order =
4001  (this)->p_order();
4002  unsigned neigh_p_order =
4004  (neigh_pt->object_pt())->p_order();
4005 
4006  // Check for p-type non-conformity
4007  if(neigh_p_order==my_p_order)
4008  {
4009  // At a conforming interface
4010  }
4011  else if(neigh_p_order<my_p_order)
4012  {
4013  // Slave at p-type non-conformity
4014  p_type_slave = true;
4015  }
4016  else
4017  {
4018  // Master at p-type non-conformity
4019  //p_type_master = true;
4020  }
4021  }
4022  // Neighbour must be smaller than me
4023  else
4024  {
4025  // Master at h-type non-conformity
4026  //h_type_master = true;
4027  }
4028  }
4029  else
4030  {
4031  // Edge is on a boundary
4032  }
4033 
4034  // Now do the mortaring
4035  //---------------------
4036  if (h_type_slave || p_type_slave)
4037  {
4038  //Compute the active coordinate index along the this side of mortar
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;
4042  else
4043  {
4044  throw OomphLibError(
4045  "Cannot transform coordinates",
4046  OOMPH_CURRENT_FUNCTION,
4047  OOMPH_EXCEPTION_LOCATION);
4048  }
4049 
4050  // Get pointer to neighbouring master element (in p-refineable form)
4052  neigh_obj_pt = dynamic_cast<PRefineableQElement<2,INITIAL_NNODE_1D>*>
4053  (neigh_pt->object_pt());
4054 
4055  // Create vector of master and slave nodes
4056  //----------------------------------------
4057  Vector<Node*> master_node_pt, slave_node_pt;
4058  Vector<unsigned> master_node_number, slave_node_number;
4059  Vector<Vector<double> > slave_node_s_fraction;
4060 
4061  //Number of nodes in one dimension
4062  const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
4063  const unsigned neigh_n_p = neigh_obj_pt->ninterpolating_node_1d(value_id);
4064 
4065  //Test for the periodic node case
4066  //Are we crossing a periodic boundary
4067  bool is_periodic = false;
4068  if(in_neighbouring_tree)
4069  {is_periodic = this->tree_pt()->root_pt()->is_neighbour_periodic(my_edge);}
4070 
4071  //If it is periodic we actually need to get the node in
4072  //the neighbour of the neighbour (which will be a parent of
4073  //the present element) so that the "fixed" coordinate is
4074  //correctly calculated.
4075  //The idea is to replace the neigh_pt and associated data
4076  //with those of the neighbour of the neighbour
4077  if(is_periodic)
4078  {
4079  throw OomphLibError(
4080  "Cannot do mortaring with periodic hanging nodes yet! (Fix me!)",
4081  "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4082  OOMPH_EXCEPTION_LOCATION);
4083  } //End of special treatment for periodic hanging nodes
4084 
4085  //Storage for pointers to the nodes and their numbers along the master edge
4086  unsigned neighbour_node_number=0;
4087  Node* neighbour_node_pt=0;
4088 
4089  // Loop over nodes along the edge
4090  for(unsigned i0=0; i0<neigh_n_p; i0++)
4091  {
4092  // Find the neighbour's node
4093  switch(neigh_edge)
4094  {
4095  case N:
4096  neighbour_node_number = i0 + neigh_n_p*(neigh_n_p-1);
4097  neighbour_node_pt
4098  = neigh_obj_pt->interpolating_node_pt(neighbour_node_number,value_id);
4099  break;
4100 
4101  case S:
4102  neighbour_node_number = i0;
4103  neighbour_node_pt
4104  = neigh_obj_pt->interpolating_node_pt(neighbour_node_number,value_id);
4105  break;
4106 
4107  case E:
4108  neighbour_node_number = neigh_n_p-1 + neigh_n_p*i0;
4109  neighbour_node_pt
4110  = neigh_obj_pt->interpolating_node_pt(neighbour_node_number,value_id);
4111  break;
4112 
4113  case W:
4114  neighbour_node_number = neigh_n_p*i0;
4115  neighbour_node_pt
4116  = neigh_obj_pt->interpolating_node_pt(neighbour_node_number,value_id);
4117  break;
4118 
4119  default:
4120  throw OomphLibError("my_edge not N, S, W, E\n",
4121  OOMPH_CURRENT_FUNCTION,
4122  OOMPH_EXCEPTION_LOCATION);
4123  }
4124 
4125  // Set node as master node
4126  master_node_number.push_back(neighbour_node_number);
4127  master_node_pt.push_back(neighbour_node_pt);
4128  }
4129 
4130  //Storage for pointers to the local nodes and their numbers along my edge
4131  unsigned local_node_number=0;
4132  Node* local_node_pt=0;
4133 
4134  // Loop over the nodes along my edge
4135  for(unsigned i0=0; i0<my_n_p; i0++)
4136  {
4137  //Storage for the fractional position of the node in the element
4138  Vector<double> s_fraction(2);
4139 
4140  // Find the local node and the fractional position of the node
4141  // which depends on the edge, of course
4142  switch(my_edge)
4143  {
4144  case N:
4145  s_fraction[0] =
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);
4150  break;
4151 
4152  case S:
4153  s_fraction[0] =
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);
4158  break;
4159 
4160  case E:
4161  s_fraction[0] = 1.0;
4162  s_fraction[1] =
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);
4166  break;
4167 
4168  case W:
4169  s_fraction[0] = 0.0;
4170  s_fraction[1] =
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);
4174  break;
4175 
4176  default:
4177  throw OomphLibError("my_edge not N, S, W, E\n",
4178  OOMPH_CURRENT_FUNCTION,
4179  OOMPH_EXCEPTION_LOCATION);
4180  }
4181 
4182  // Add node to vector of slave element nodes
4183  slave_node_number.push_back(local_node_number);
4184  slave_node_pt.push_back(local_node_pt);
4185 
4186  // Store node's local fraction
4187  slave_node_s_fraction.push_back(s_fraction);
4188  }
4189 
4190  // Store the number of slave and master nodes for use later
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;
4195 
4196  // Storage for master shapes
4197  Shape master_shapes(neigh_obj_pt->ninterpolating_node(value_id));
4198 
4199  // Get master and slave nodal positions
4200  //-------------------------------------
4201  Vector<double> slave_nodal_position;
4202  Vector<double> slave_weight(slave_element_nnode_1d);
4203  Orthpoly::gll_nodes(slave_element_nnode_1d,
4204  slave_nodal_position,slave_weight);
4205  Vector<double> master_nodal_position;
4206  Vector<double> master_weight(master_element_nnode_1d);
4207  Orthpoly::gll_nodes(master_element_nnode_1d,
4208  master_nodal_position,master_weight);
4209 
4210  // Apply the vertex matching condition
4211  //------------------------------------
4212  // Vertiex matching is ensured automatically in cases where there is a node
4213  // at each end of the mortar that is shared between the master and slave
4214  // elements. Where this is not the case, the vertex matching condition must
4215  // be enforced by constraining the value of the unknown at the node on the
4216  // slave side to be the same as the value at that location in the master.
4217 
4218  // Store positions of the mortar vertex/non-vertex nodes in the slave element
4219  const unsigned n_mortar_vertices = 2;
4220  Vector<unsigned> vertex_pos(n_mortar_vertices);
4221  vertex_pos[0] = 0;
4222  vertex_pos[1] = this->ninterpolating_node_1d(value_id)-1;
4223  Vector<unsigned> non_vertex_pos(my_n_p-n_mortar_vertices);
4224  for(unsigned i=0; i<my_n_p-n_mortar_vertices; i++)
4225  {non_vertex_pos[i] = i+1;}
4226 
4227  // Check if the mortar vertices are shared
4228  for (unsigned v=0; v<n_mortar_vertices; v++)
4229  {
4230  // Search master node storage for the node
4231  bool node_is_shared=false;
4232  for(unsigned i=0; i<master_node_pt.size(); i++)
4233  {
4234  if(slave_node_pt[vertex_pos[v]]==master_node_pt[i])
4235  {
4236  node_is_shared=true;
4237  break;
4238  }
4239  }
4240 
4241  // If the node is not shared then we must constrain its value by setting
4242  // up a hanging scheme
4243  if(!node_is_shared)
4244  {
4245  // Calculate weights. These are just the master shapes evaluated at
4246  // this slave node's position
4247 
4248  // Work out this node's location in the master
4249  Vector<double> s_in_neigh(2);
4250  for(unsigned i=0;i<2;i++)
4251  {
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]);
4255  }
4256 
4257  // Get master shapes at slave nodal positions
4258  neigh_obj_pt->interpolating_basis(s_in_neigh,master_shapes,value_id);
4259 
4260  // Remove any existing hanging node info
4261  // (This may be a bit wasteful, but guarantees correctness)
4262  slave_node_pt[vertex_pos[v]]->set_nonhanging();
4263 
4264  // Set up hanging scheme for this node
4265  HangInfo* explicit_hang_pt = new HangInfo(n_master_nodes);
4266  for(unsigned m=0; m<n_master_nodes; m++)
4267  {
4268  explicit_hang_pt->set_master_node_pt(m,master_node_pt[m],master_shapes[master_node_number[m]]);
4269  }
4270 
4271  // Make node hang
4272  slave_node_pt[vertex_pos[v]]->set_hanging_pt(explicit_hang_pt,-1);
4273  }
4274  }
4275 
4276  // Check that there are actually slave nodes for which we still need to
4277  // construct a hanging scheme. If not then there is nothing more to do.
4278  if(n_slave_nodes-n_mortar_vertices > 0)
4279  {
4280 
4281  // Assemble mass matrix for mortar
4282  //--------------------------------
4283  Vector<double> psi(n_slave_nodes-n_mortar_vertices);
4284  Vector<double> diag_M(n_slave_nodes-n_mortar_vertices);
4285  Vector<Vector<double> > shared_node_M(n_mortar_vertices);
4286  for (unsigned i=0; i<shared_node_M.size(); i++)
4287  {shared_node_M[i].resize(n_slave_nodes-n_mortar_vertices);}
4288 
4289  // Fill in part corresponding to slave nodal positions (unknown)
4290  for (unsigned i=0; i<n_slave_nodes-n_mortar_vertices; i++)
4291  {
4292  // Use L'Hosptal's rule:
4293  psi[i] = pow(-1.0,int((slave_element_nnode_1d-1)-i-1))
4294  *-Orthpoly::ddlegendre(slave_element_nnode_1d-1,
4295  slave_nodal_position[non_vertex_pos[i]]);
4296  // Put in contribution
4297  diag_M[i] = psi[i]*slave_weight[non_vertex_pos[i]];
4298  }
4299 
4300  // Fill in part corresponding to slave element vertices (known)
4301  for(unsigned v=0; v<shared_node_M.size(); v++)
4302  {
4303  for (unsigned i=0; i<n_slave_nodes-n_mortar_vertices; i++)
4304  {
4305  // Check if denominator is zero
4306  if (std::fabs(slave_nodal_position[non_vertex_pos[i]]
4307  - slave_nodal_position[vertex_pos[v]]) >= 1.0e-8 )
4308  {
4309  // We're ok
4310  psi[i] = pow(-1.0,int((slave_element_nnode_1d-1)-i-1))
4311  * Orthpoly::dlegendre(slave_element_nnode_1d-1,
4312  slave_nodal_position[vertex_pos[v]])
4313  / (slave_nodal_position[non_vertex_pos[i]]
4314  - slave_nodal_position[vertex_pos[v]]);
4315  }
4316  // Check if numerator is zero
4317  else if (std::fabs(Orthpoly::dlegendre(slave_element_nnode_1d-1,
4318  slave_nodal_position[vertex_pos[v]]))
4319  < 1.0e-8)
4320  {
4321  // We can use l'hopital's rule
4322  psi[i] = pow(-1.0,int((slave_element_nnode_1d-1)-i-1))
4323  *-Orthpoly::ddlegendre(slave_element_nnode_1d-1,
4324  slave_nodal_position[non_vertex_pos[i]]);
4325  }
4326  else
4327  {
4328  // We can't use l'hopital's rule
4329  throw
4330  OomphLibError(
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);
4334  }
4335  // Put in contribution
4336  shared_node_M[v][i] = psi[i]*slave_weight[vertex_pos[v]];
4337  }
4338  }
4339 
4340  // Assemble local projection matrix for mortar
4341  //--------------------------------------------
4342 
4343  // Have only one local projection matrix because there is just one master
4344  Vector<Vector<double> > P(n_slave_nodes-n_mortar_vertices);
4345  for(unsigned i=0; i<P.size(); i++) {P[i].resize(n_master_nodes,0.0);}
4346 
4347  //Storage for local coordinate
4348  Vector<double> s(2);
4349 
4350  // Sum contributions from master element shapes (quadrature).
4351  // The order of quadrature must be high enough to evaluate a polynomial
4352  // of order N_s + N_m - 3 exactly, where N_s = n_slave_nodes, N_m =
4353  // n_master_nodes.
4354  // (Use pointers for the quadrature knots and weights so that
4355  // data is not unnecessarily copied)
4356  //unsigned quadrature_order =
4357  // std::max(slave_element_nnode_1d,master_element_nnode_1d);
4358  Vector<double> *quadrature_knot, *quadrature_weight;
4359  if (slave_element_nnode_1d >= master_element_nnode_1d)
4360  {
4361  // Use the same quadrature order as the slave element (me)
4362  quadrature_knot = &slave_nodal_position;
4363  quadrature_weight = &slave_weight;
4364  }
4365  else
4366  {
4367  // Use the same quadrature order as the master element (neighbour)
4368  quadrature_knot = &master_nodal_position;
4369  quadrature_weight = &master_weight;
4370  }
4371 
4372  // Quadrature loop
4373  for (unsigned q=0; q<(*quadrature_weight).size(); q++)
4374  {
4375  // Evaluate mortar test functions at the quadrature knots in the slave
4376  //s[active_coord_index] = (*quadrature_knot)[q];
4377  Vector<double> s_on_mortar(1);
4378  s_on_mortar[0] = (*quadrature_knot)[q];
4379 
4380  // Get psi
4381  for(unsigned k=0; k<n_slave_nodes-n_mortar_vertices; k++)
4382  {
4383  // Check if denominator is zero
4384  if (std::fabs(slave_nodal_position[non_vertex_pos[k]]-s_on_mortar[0]) >= 1.0e-08)
4385  {
4386  // We're ok
4387  psi[k] = pow(-1.0,int((slave_element_nnode_1d-1)-k-1))
4388  * Orthpoly::dlegendre(slave_element_nnode_1d-1,s_on_mortar[0])
4389  / (slave_nodal_position[non_vertex_pos[k]]-s_on_mortar[0]);
4390  }
4391  // Check if numerator is zero
4392  else if (std::fabs(Orthpoly::dlegendre(slave_element_nnode_1d-1,s_on_mortar[0]))
4393  < 1.0e-8)
4394  {
4395  // We can use l'Hopital's rule
4396  psi[k] = pow(-1.0,int((slave_element_nnode_1d-1)-k-1))
4397  * -Orthpoly::ddlegendre(slave_element_nnode_1d-1,s_on_mortar[0]);
4398  }
4399  else
4400  {
4401  // We can't use l'hopital's rule
4402  throw
4403  OomphLibError(
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);
4407  }
4408  }
4409 
4410  // Convert coordinate on mortar to local fraction in slave element
4411  Vector<double> s_fraction(2);
4412  for(unsigned i=0; i<2; i++)
4413  {
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];
4417  }
4418 
4419  // Project active coordinate into master element
4420  Vector<double> s_in_neigh(2);
4421  for(unsigned i=0;i<2;i++)
4422  {
4423  s_in_neigh[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
4424  (s_hi_neigh[i] - s_lo_neigh[i]);
4425  }
4426 
4427  // Evaluate master shapes at projections of local quadrature knots
4428  neigh_obj_pt->interpolating_basis(s_in_neigh,master_shapes,value_id);
4429 
4430  // Populate local projection matrix
4431  for(unsigned i=0; i<n_slave_nodes-n_mortar_vertices; i++)
4432  {
4433  for(unsigned j=0; j<n_master_nodes; j++)
4434  {
4435  P[i][j] += master_shapes[master_node_number[j]]*psi[i]*(*quadrature_weight)[q];
4436  }
4437  }
4438  }
4439 
4440  // Assemble global projection matrices for mortar
4441  //-----------------------------------------------
4442  // Need to subtract contributions from the "known unknowns" corresponding
4443  // to the nodes at the vertices of the mortar
4444 
4445  // Assemble contributions from mortar vertex nodes
4446  for(unsigned v=0; v<n_mortar_vertices; v++)
4447  {
4448  // Convert coordinate on mortar to local fraction in slave element
4449  Vector<double> s_fraction(2);
4450  for(unsigned i=0; i<2; i++)
4451  {
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];
4455  }
4456 
4457  // Project active coordinate into master element
4458  Vector<double> s_in_neigh(2);
4459  for(unsigned i=0;i<2;i++)
4460  {
4461  s_in_neigh[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
4462  (s_hi_neigh[i] - s_lo_neigh[i]);
4463  }
4464 
4465  // Get master shapes at slave nodal positions
4466  neigh_obj_pt->interpolating_basis(s_in_neigh,master_shapes,value_id);
4467 
4468  for(unsigned i=0; i<n_slave_nodes-n_mortar_vertices; i++)
4469  {
4470  for(unsigned k=0; k<n_master_nodes; k++)
4471  {
4472  P[i][k] -= master_shapes[master_node_number[k]]*shared_node_M[v][i];
4473  }
4474  }
4475  }
4476 
4477  // Solve mortar system
4478  //--------------------
4479  for(unsigned i=0; i<n_slave_nodes-n_mortar_vertices; i++)
4480  {
4481  for(unsigned j=0; j<n_master_nodes; j++)
4482  {
4483  P[i][j]/=diag_M[i];
4484  }
4485  }
4486 
4487  // Create structures to hold the hanging info
4488  //-------------------------------------------
4489  Vector<HangInfo*> hang_info_pt(n_slave_nodes);
4490  for (unsigned i=0; i<n_slave_nodes-n_mortar_vertices; i++)
4491  {
4492  hang_info_pt[i] = new HangInfo(n_master_nodes);
4493  }
4494 
4495  // Copy information to hanging nodes
4496  //----------------------------------
4497  for(unsigned i=0; i<n_slave_nodes-n_mortar_vertices; i++)
4498  {
4499  for(unsigned j=0; j<n_master_nodes; j++)
4500  {
4501  hang_info_pt[i]->set_master_node_pt(j,master_node_pt[j],P[i][j]);
4502  }
4503  }
4504 
4505  // Set pointers to hanging info
4506  //-----------------------------
4507  for (unsigned i=0; i<n_slave_nodes-n_mortar_vertices; i++)
4508  {
4509  // Check that the node shoule actually hang.
4510  // That is, if the polynomial orders of the elements at a p-type
4511  // non-conormity are both odd then the middle node on the edge is
4512  // shared but a hanging scheme has been constructed for it.
4513  bool node_is_really_shared = false;
4514  for (unsigned m=0; m<hang_info_pt[i]->nmaster(); m++)
4515  {
4516  // We can simply check if the hanging scheme lists itself as a master
4517  if (hang_info_pt[i]->master_node_pt(m)==slave_node_pt[non_vertex_pos[i]])
4518  {
4519  node_is_really_shared = true;
4520 
4521 #ifdef PARANOID
4522  // Also check the corresponding weight: it should be 1
4523  if(std::fabs(hang_info_pt[i]->master_weight(m)-1.0) > 1.0e-06)
4524  {
4525  throw
4526  OomphLibError(
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);
4530  }
4531 #endif
4532  }
4533  }
4534 
4535  // Now we can make the node hang, if it isn't a shared node
4536  if(!node_is_really_shared)
4537  {
4538  slave_node_pt[non_vertex_pos[i]]->set_hanging_pt(hang_info_pt[i],-1);
4539  }
4540  }
4541 
4542  } // End of case where there are still slave nodes
4543 
4544 #ifdef PARANOID
4545  // Check all slave nodes, hanging or otherwise
4546  for (unsigned i=0; i<n_slave_nodes; i++)
4547  {
4548  // Check that weights sum to 1 for those that hang
4549  if (slave_node_pt[i]->is_hanging())
4550  {
4551  // Check that weights sum to 1
4552  double weight_sum = 0.0;
4553  for (unsigned m=0; m<slave_node_pt[i]->hanging_pt()->nmaster(); m++)
4554  {
4555  weight_sum += slave_node_pt[i]->hanging_pt()->master_weight(m);
4556  }
4557 
4558  // Warn if not
4559  if (std::fabs(weight_sum-1.0) > 1.0e-08)
4560  {
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);
4566  }
4567  }
4568  else
4569  {
4570  // Check that this node is shared with the master element if it
4571  // isn't hanging
4572  bool is_master = false;
4573  for (unsigned n=0; n<n_master_nodes; n++)
4574  {
4575  if (slave_node_pt[i]==master_node_pt[n])
4576  {
4577  // Node is a master
4578  is_master = true;
4579  break;
4580  }
4581  }
4582 
4583  if (!is_master)
4584  {
4585  // Throw error
4586  std::ostringstream error_string;
4587  error_string
4588  << "This node in the slave element is neither" << std::endl
4589  << "hanging or shared with a master element." << std::endl;
4590 
4591  throw OomphLibError(
4592  error_string.str(),
4593  "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4594  OOMPH_EXCEPTION_LOCATION);
4595  }
4596  }
4597  }
4598 #endif
4599 
4600  // Finally, Loop over all slave nodes and fine-tune their positions
4601  //-----------------------------------------------------------------
4602  // Here we simply set the node's positions to be consistent
4603  // with the hanging scheme. This is not strictly necessary
4604  // because it is done in the mesh adaptation before the node
4605  // becomes non-hanging later on. We make no attempt to ensure
4606  // (strong) continuity in the position across the mortar.
4607  for(unsigned i=0; i<n_slave_nodes; i++)
4608  {
4609  // Only fine-tune hanging nodes
4610  if(slave_node_pt[i]->is_hanging())
4611  {
4612  //If we are doing the position, then
4613  if(value_id==-1)
4614  {
4615  // Get the local coordinate of this slave node
4616  Vector<double> s_local(2);
4617  this->local_coordinate_of_node(slave_node_number[i],s_local);
4618 
4619  // Get the position from interpolation in this element via
4620  // the hanging scheme
4621  Vector<double> x_in_neighb(2);
4622  this->interpolated_x(s_local,x_in_neighb);
4623 
4624  // Fine adjust the coordinates (macro map will pick up boundary
4625  // accurately but will lead to different element edges)
4626  slave_node_pt[i]->x(0)=x_in_neighb[0];
4627  slave_node_pt[i]->x(1)=x_in_neighb[1];
4628  }
4629  }
4630  }
4631  } //End of case where this interface is to be mortared
4632 
4633 }
4634 
4635 ////////////////////////////////////////////////////////////////
4636 // 3D PRefineableQElements
4637 ////////////////////////////////////////////////////////////////
4638 
4639 /// Get local coordinates of node j in the element; vector sets its own size
4640 template<unsigned INITIAL_NNODE_1D>
4642 local_coordinate_of_node(const unsigned& n, Vector<double>& s) const
4643  {
4644  s.resize(3);
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));
4649 
4650  switch(Nnode_1d)
4651  {
4652  case 2:
4657  break;
4658  case 3:
4663  break;
4664  case 4:
4669  break;
4670  case 5:
4675  break;
4676  case 6:
4681  break;
4682  case 7:
4687  break;
4688  default:
4689  std::ostringstream error_message;
4690  error_message <<"\nERROR: Exceeded maximum polynomial order for";
4691  error_message <<"\n shape functions.\n";
4692  throw OomphLibError(error_message.str(),
4693  OOMPH_CURRENT_FUNCTION,
4694  OOMPH_EXCEPTION_LOCATION);
4695  break;
4696  }
4697  }
4698 
4699 /// Get the local fractino of node j in the element
4700 template<unsigned INITIAL_NNODE_1D>
4702 local_fraction_of_node(const unsigned &n, Vector<double> &s_fraction)
4703  {
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);
4708  }
4709 
4710 template<unsigned INITIAL_NNODE_1D>
4712 local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
4713  {
4714  switch(this->nnode_1d())
4715  {
4716  case 2:
4718  return 0.5*(OneDimensionalLegendreShape<2>::nodal_position(n1d) + 1.0);
4719  case 3:
4721  return 0.5*(OneDimensionalLegendreShape<3>::nodal_position(n1d) + 1.0);
4722  case 4:
4724  return 0.5*(OneDimensionalLegendreShape<4>::nodal_position(n1d) + 1.0);
4725  case 5:
4727  return 0.5*(OneDimensionalLegendreShape<5>::nodal_position(n1d) + 1.0);
4728  case 6:
4730  return 0.5*(OneDimensionalLegendreShape<6>::nodal_position(n1d) + 1.0);
4731  case 7:
4733  return 0.5*(OneDimensionalLegendreShape<7>::nodal_position(n1d) + 1.0);
4734  default:
4735  std::ostringstream error_message;
4736  error_message <<"\nERROR: Exceeded maximum polynomial order for";
4737  error_message <<"\n shape functions.\n";
4738  throw OomphLibError(error_message.str(),
4739  OOMPH_CURRENT_FUNCTION,
4740  OOMPH_EXCEPTION_LOCATION);
4741  return 0.0;
4742  }
4743  }
4744 
4745 //==================================================================
4746 /// Return the node at the specified local coordinate
4747 //==================================================================
4748 template<unsigned INITIAL_NNODE_1D>
4751 {
4752  unsigned Nnode_1d = this->nnode_1d();
4753  //Load the tolerance into a local variable
4755  //There are two possible indices.
4756  Vector<int> index(3);
4757 
4758  // Loop over indices
4759  for (unsigned i=0; i<3; i++)
4760  {
4761  // Determine the index
4762  // -------------------
4763 
4764  bool is_found=false;
4765 
4766  // If we are at the lower limit, the index is zero
4767  if(std::fabs(s[i] + 1.0) < tol)
4768  {
4769  index[i] = 0;
4770  is_found=true;
4771  }
4772  // If we are at the upper limit, the index is the number of nodes minus 1
4773  else if(std::fabs(s[i] - 1.0) < tol)
4774  {
4775  index[i] = Nnode_1d-1;
4776  is_found=true;
4777  }
4778  // Otherwise, we have to calculate the index in general
4779  else
4780  {
4781  // Compute Gauss-Lobatto-Legendre node positions
4782  Vector<double> z;
4783  Orthpoly::gll_nodes(Nnode_1d, z);
4784  // Loop over possible internal nodal positions
4785  for (unsigned n=1; n<Nnode_1d-1; n++)
4786  {
4787  if (std::fabs(z[n] - s[i]) < tol)
4788  {
4789  index[i] = n;
4790  is_found=true;
4791  break;
4792  }
4793  }
4794  }
4795 
4796  if (!is_found)
4797  {
4798  // No matching nodes
4799  return 0;
4800  }
4801  }
4802  // If we've got here we have a node, so let's return a pointer to it
4803  return this->node_pt(index[0] + Nnode_1d*index[1] + Nnode_1d*Nnode_1d*index[2]);
4804 }
4805 
4806 //===================================================================
4807 /// If a neighbouring element has already created a node at
4808 /// a position corresponding to the local fractional position within the
4809 /// present element, s_fraction, return
4810 /// a pointer to that node. If not, return NULL (0).
4811 //===================================================================
4812 template<unsigned INITIAL_NNODE_1D>
4815 {
4816  using namespace OcTreeNames;
4817 
4818  //Calculate the faces/edges on which the node lies
4819  Vector<int> faces;
4820  Vector<int> edges;
4821 
4822  if(s_fraction[0]==0.0)
4823  {
4824  faces.push_back(L);
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);}
4829  }
4830 
4831  if(s_fraction[0]==1.0)
4832  {
4833  faces.push_back(R);
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);}
4838  }
4839 
4840  if(s_fraction[1]==0.0)
4841  {
4842  faces.push_back(D);
4843  if (s_fraction[2]==0.0) {edges.push_back(DB);}
4844  if (s_fraction[2]==1.0) {edges.push_back(DF);}
4845  }
4846 
4847  if(s_fraction[1]==1.0)
4848  {
4849  faces.push_back(U);
4850  if (s_fraction[2]==0.0) {edges.push_back(UB);}
4851  if (s_fraction[2]==1.0) {edges.push_back(UF);}
4852  }
4853 
4854  if(s_fraction[2]==0.0)
4855  {
4856  faces.push_back(B);
4857  }
4858 
4859  if(s_fraction[2]==1.0)
4860  {
4861  faces.push_back(F);
4862  }
4863 
4864  //Find the number of faces
4865  unsigned n_face = faces.size();
4866 
4867  //Find the number of edges
4868  unsigned n_edge = edges.size();
4869 
4870  Vector<unsigned> translate_s(3);
4871  Vector<double> s_lo_neigh(3);
4872  Vector<double> s_hi_neigh(3);
4873  Vector<double> s(3);
4874 
4875  int neigh_face, diff_level;
4876 
4877  //Loop over the faces on which the node lies
4878  //------------------------------------------
4879  for(unsigned j=0;j<n_face;j++)
4880  {
4881  // Find pointer to neighbouring element along face
4882  OcTree* neigh_pt;
4883  neigh_pt = octree_pt()->
4884  gteq_face_neighbour(faces[j],translate_s,s_lo_neigh,s_hi_neigh,neigh_face,
4885  diff_level);
4886 
4887  // Neighbour exists
4888  if(neigh_pt!=0)
4889  {
4890  // Have any of its vertex nodes been created yet?
4891  // (Must look in incomplete neighbours because after the
4892  // pre-build they may contain pointers to the required nodes. e.g.
4893  // h-refinement of neighbouring linear and quadratic elements)
4894  bool a_vertex_node_is_built = false;
4895  QElement<3,INITIAL_NNODE_1D>* neigh_obj_pt =
4896  dynamic_cast<QElement<3,INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
4897  if(neigh_obj_pt==0)
4898  {
4899  throw
4900  OomphLibError("Not a quad element!",
4901  "PRefineableQElement<3,INITIAL_NNODE_1D>::node_created_by_neighbour()",
4902  OOMPH_EXCEPTION_LOCATION);
4903  }
4904  for(unsigned vnode=0; vnode<neigh_obj_pt->nvertex_node(); vnode++)
4905  {
4906  if(neigh_obj_pt->vertex_node_pt(vnode)!=0)
4907  a_vertex_node_is_built = true;
4908  break;
4909  }
4910  if(a_vertex_node_is_built)
4911  {
4912  //We now need to translate the nodal location, defined in terms
4913  //of the fractional coordinates of the present element into
4914  //those of its neighbour. For this we use the information returned
4915  //to use from the octree function.
4916 
4917  //Calculate the local coordinate in the neighbour
4918  //Note that we need to use the translation scheme in case
4919  //the local coordinates are swapped in the neighbour.
4920  for(unsigned i=0;i<3;i++)
4921  {
4922  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
4923  (s_hi_neigh[i] - s_lo_neigh[i]);
4924  }
4925 
4926  //Find the node in the neighbour
4927  Node* neighbour_node_pt =
4928  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
4929 
4930  //If there is a node, return it
4931  if(neighbour_node_pt!=0)
4932  {
4933  return neighbour_node_pt;
4934  }
4935  }
4936  }
4937  } //End of loop over faces
4938 
4939 
4940  //Loop over the edges on which the node lies
4941  //------------------------------------------
4942  for(unsigned j=0;j<n_edge;j++)
4943  {
4944 
4945  // Even if we restrict ourselves to true edge neighbours (i.e.
4946  // elements that are not also face neighbours) there may be multiple
4947  // edge neighbours across the edges between multiple root octrees.
4948  // When making the first call to OcTree::gteq_true_edge_neighbour(...)
4949  // we simply return the first one of these multiple edge neighbours
4950  // (if there are any at all, of course) and also return the total number
4951  // of true edge neighbours. If the node in question already exists
4952  // on the first edge neighbour we're done. If it doesn't it may exist
4953  // on other edge neighbours so we repeat the process over all
4954  // other edge neighbours (bailing out if a node is found, of course).
4955 
4956  // Initially return the zero-th true edge neighbour
4957  unsigned i_root_edge_neighbour=0;
4958 
4959  // Initialise the total number of true edge neighbours
4960  unsigned nroot_edge_neighbour=0;
4961 
4962  // Keep searching until we've found the node or until we've checked
4963  // all available edge neighbours
4964  bool keep_searching=true;
4965  while (keep_searching)
4966  {
4967 
4968  // Find pointer to neighbouring element along edge
4969  OcTree* neigh_pt;
4970  neigh_pt = octree_pt()->
4971  gteq_true_edge_neighbour(edges[j],
4972  i_root_edge_neighbour,
4973  nroot_edge_neighbour,
4974  translate_s,
4975  s_lo_neigh,s_hi_neigh,neigh_face,
4976  diff_level);
4977 
4978  // Neighbour exists
4979  if(neigh_pt!=0)
4980  {
4981  // Have any of its vertex nodes been created yet?
4982  // (Must look in incomplete neighbours because after the
4983  // pre-build they may contain pointers to the required nodes. e.g.
4984  // h-refinement of neighbouring linear and quadratic elements)
4985  bool a_vertex_node_is_built = false;
4986  QElement<3,INITIAL_NNODE_1D>* neigh_obj_pt =
4987  dynamic_cast<QElement<3,INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
4988  if(neigh_obj_pt==0)
4989  {
4990  throw
4991  OomphLibError("Not a quad element!",
4992  "PRefineableQElement<3,INITIAL_NNODE_1D>::node_created_by_neighbour()",
4993  OOMPH_EXCEPTION_LOCATION);
4994  }
4995  for(unsigned vnode=0; vnode<neigh_obj_pt->nvertex_node(); vnode++)
4996  {
4997  if(neigh_obj_pt->vertex_node_pt(vnode)!=0)
4998  a_vertex_node_is_built = true;
4999  break;
5000  }
5001  if(a_vertex_node_is_built)
5002  {
5003  //We now need to translate the nodal location, defined in terms
5004  //of the fractional coordinates of the present element into
5005  //those of its neighbour. For this we use the information returned
5006  //to use from the octree function.
5007 
5008  //Calculate the local coordinate in the neighbour
5009  //Note that we need to use the translation scheme in case
5010  //the local coordinates are swapped in the neighbour.
5011  for(unsigned i=0;i<3;i++)
5012  {
5013  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
5014  (s_hi_neigh[i] - s_lo_neigh[i]);
5015  }
5016 
5017  //Find the node in the neighbour
5018  Node* neighbour_node_pt =
5019  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
5020 
5021  //If there is a node, return it
5022  if(neighbour_node_pt!=0)
5023  {
5024  return neighbour_node_pt;
5025  }
5026  }
5027  }
5028 
5029  // Keep searching, but only if there are further edge neighbours
5030  // Try next root edge neighbour
5031  i_root_edge_neighbour++;
5032 
5033  // Have we exhausted the search
5034  if (i_root_edge_neighbour>=nroot_edge_neighbour)
5035  {
5036  keep_searching=false;
5037  }
5038 
5039  } // End of while keep searching over all true edge neighbours
5040 
5041  } //End of loop over edges
5042 
5043  //Node not found, return null
5044  return 0;
5045 }
5046 
5047 //===================================================================
5048 /// If a neighbouring element's son has already created a node at
5049 /// a position corresponding to the local fractional position within the
5050 /// present element, s_fraction, return
5051 /// a pointer to that node. If not, return NULL (0). If the node is
5052 /// periodic the flag is_periodic will be true
5053 //===================================================================
5054 template<unsigned INITIAL_NNODE_1D>
5057 {
5058  using namespace OcTreeNames;
5059 
5060  //Calculate the faces/edges on which the node lies
5061  Vector<int> faces;
5062  Vector<int> edges;
5063 
5064  if(s_fraction[0]==0.0)
5065  {
5066  faces.push_back(L);
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);}
5071  }
5072 
5073  if(s_fraction[0]==1.0)
5074  {
5075  faces.push_back(R);
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);}
5080  }
5081 
5082  if(s_fraction[1]==0.0)
5083  {
5084  faces.push_back(D);
5085  if (s_fraction[2]==0.0) {edges.push_back(DB);}
5086  if (s_fraction[2]==1.0) {edges.push_back(DF);}
5087  }
5088 
5089  if(s_fraction[1]==1.0)
5090  {
5091  faces.push_back(U);
5092  if (s_fraction[2]==0.0) {edges.push_back(UB);}
5093  if (s_fraction[2]==1.0) {edges.push_back(UF);}
5094  }
5095 
5096  if(s_fraction[2]==0.0)
5097  {
5098  faces.push_back(B);
5099  }
5100 
5101  if(s_fraction[2]==1.0)
5102  {
5103  faces.push_back(F);
5104  }
5105 
5106  //Find the number of faces
5107  unsigned n_face = faces.size();
5108 
5109  //Find the number of edges
5110  unsigned n_edge = edges.size();
5111 
5112  Vector<unsigned> translate_s(3);
5113  Vector<double> s_lo_neigh(3);
5114  Vector<double> s_hi_neigh(3);
5115  Vector<double> s(3);
5116 
5117  int neigh_face, diff_level;
5118 
5119  //Loop over the faces on which the node lies
5120  //------------------------------------------
5121  for(unsigned j=0;j<n_face;j++)
5122  {
5123  // Find pointer to neighbouring element along face
5124  OcTree* neigh_pt;
5125  neigh_pt = octree_pt()->
5126  gteq_face_neighbour(faces[j],translate_s,s_lo_neigh,s_hi_neigh,neigh_face,
5127  diff_level);
5128 
5129  // Neighbour exists
5130  if(neigh_pt!=0)
5131  {
5132  // Have its nodes been created yet?
5133  // (Must look in sons of unfinished neighbours too!!!)
5134  if(true)
5135  {
5136  //We now need to translate the nodal location, defined in terms
5137  //of the fractional coordinates of the present element into
5138  //those of its neighbour. For this we use the information returned
5139  //to use from the octree function.
5140 
5141  //Calculate the local coordinate in the neighbour
5142  //Note that we need to use the translation scheme in case
5143  //the local coordinates are swapped in the neighbour.
5144  for(unsigned i=0;i<3;i++)
5145  {
5146  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
5147  (s_hi_neigh[i] - s_lo_neigh[i]);
5148  }
5149 
5150  // Check if the element has sons
5151  if(neigh_pt->nsons()!=0)
5152  {
5153  //First, find the son element in which the node should live
5154 
5155  //Find coordinates in the sons
5156  Vector<double> s_in_son(3);
5157  int son=-10;
5158  //On the left
5159  if(s[0] < 0.0)
5160  {
5161  //On the bottom
5162  if(s[1] < 0.0)
5163  {
5164  //On the back
5165  if(s[2] < 0.0)
5166  {
5167  //It's the LDB son
5168  son = OcTreeNames::LDB;
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];
5172  }
5173  //On the front
5174  else
5175  {
5176  //It's the LDF son
5177  son = OcTreeNames::LDF;
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];
5181  }
5182  }
5183  //On the top
5184  else
5185  {
5186  //On the back
5187  if(s[2] < 0.0)
5188  {
5189  //It's the LUB son
5190  son = OcTreeNames::LUB;
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];
5194  }
5195  //On the front
5196  else
5197  {
5198  //It's the LUF son
5199  son = OcTreeNames::LUF;
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];
5203  }
5204  }
5205  }
5206  //On the right
5207  else
5208  {
5209  //On the bottom
5210  if(s[1] < 0.0)
5211  {
5212  //On the back
5213  if(s[2] < 0.0)
5214  {
5215  //It's the RDB son
5216  son = OcTreeNames::RDB;
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];
5220  }
5221  //On the front
5222  else
5223  {
5224  //It's the RDF son
5225  son = OcTreeNames::RDF;
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];
5229  }
5230  }
5231  //On the top
5232  else
5233  {
5234  //On the back
5235  if(s[2] < 0.0)
5236  {
5237  //It's the RUB son
5238  son = OcTreeNames::RUB;
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];
5242  }
5243  //On the front
5244  else
5245  {
5246  //It's the RUF son
5247  son = OcTreeNames::RUF;
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];
5251  }
5252  }
5253  }
5254 
5255  //Find the node in the neighbour's son
5256  Node* neighbour_son_node_pt =
5257  neigh_pt->son_pt(son)->object_pt()->
5258  get_node_at_local_coordinate(s_in_son);
5259 
5260  //If there is a node, return it
5261  if(neighbour_son_node_pt!=0)
5262  {
5263  //Return the pointer to the neighbouring node
5264  return neighbour_son_node_pt;
5265  }
5266  }
5267  }
5268  }
5269  } //End of loop over faces
5270 
5271 
5272  //Loop over the edges on which the node lies
5273  //------------------------------------------
5274  for(unsigned j=0;j<n_edge;j++)
5275  {
5276 
5277  // Even if we restrict ourselves to true edge neighbours (i.e.
5278  // elements that are not also face neighbours) there may be multiple
5279  // edge neighbours across the edges between multiple root octrees.
5280  // When making the first call to OcTree::gteq_true_edge_neighbour(...)
5281  // we simply return the first one of these multiple edge neighbours
5282  // (if there are any at all, of course) and also return the total number
5283  // of true edge neighbours. If the node in question already exists
5284  // on the first edge neighbour we're done. If it doesn't it may exist
5285  // on other edge neighbours so we repeat the process over all
5286  // other edge neighbours (bailing out if a node is found, of course).
5287 
5288  // Initially return the zero-th true edge neighbour
5289  unsigned i_root_edge_neighbour=0;
5290 
5291  // Initialise the total number of true edge neighbours
5292  unsigned nroot_edge_neighbour=0;
5293 
5294  // Keep searching until we've found the node or until we've checked
5295  // all available edge neighbours
5296  bool keep_searching=true;
5297  while (keep_searching)
5298  {
5299 
5300  // Find pointer to neighbouring element along edge
5301  OcTree* neigh_pt;
5302  neigh_pt = octree_pt()->
5303  gteq_true_edge_neighbour(edges[j],
5304  i_root_edge_neighbour,
5305  nroot_edge_neighbour,
5306  translate_s,
5307  s_lo_neigh,s_hi_neigh,neigh_face,
5308  diff_level);
5309 
5310  // Neighbour exists
5311  if(neigh_pt!=0)
5312  {
5313  // Have its nodes been created yet?
5314  // (Must look in sons of unfinished neighbours too!!!)
5315  if(true)
5316  {
5317  //We now need to translate the nodal location, defined in terms
5318  //of the fractional coordinates of the present element into
5319  //those of its neighbour. For this we use the information returned
5320  //to use from the octree function.
5321 
5322  //Calculate the local coordinate in the neighbour
5323  //Note that we need to use the translation scheme in case
5324  //the local coordinates are swapped in the neighbour.
5325  for(unsigned i=0;i<3;i++)
5326  {
5327  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
5328  (s_hi_neigh[i] - s_lo_neigh[i]);
5329  }
5330 
5331  // Check if the element has sons
5332  if(neigh_pt->nsons()!=0)
5333  {
5334  //First, find the son element in which the node should live
5335 
5336  //Find coordinates in the sons
5337  Vector<double> s_in_son(3);
5338  int son=-10;
5339  //On the left
5340  if(s[0] < 0.0)
5341  {
5342  //On the bottom
5343  if(s[1] < 0.0)
5344  {
5345  //On the back
5346  if(s[2] < 0.0)
5347  {
5348  //It's the LDB son
5349  son = OcTreeNames::LDB;
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];
5353  }
5354  //On the front
5355  else
5356  {
5357  //It's the LDF son
5358  son = OcTreeNames::LDF;
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];
5362  }
5363  }
5364  //On the top
5365  else
5366  {
5367  //On the back
5368  if(s[2] < 0.0)
5369  {
5370  //It's the LUB son
5371  son = OcTreeNames::LUB;
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];
5375  }
5376  //On the front
5377  else
5378  {
5379  //It's the LUF son
5380  son = OcTreeNames::LUF;
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];
5384  }
5385  }
5386  }
5387  //On the right
5388  else
5389  {
5390  //On the bottom
5391  if(s[1] < 0.0)
5392  {
5393  //On the back
5394  if(s[2] < 0.0)
5395  {
5396  //It's the RDB son
5397  son = OcTreeNames::RDB;
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];
5401  }
5402  //On the front
5403  else
5404  {
5405  //It's the RDF son
5406  son = OcTreeNames::RDF;
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];
5410  }
5411  }
5412  //On the top
5413  else
5414  {
5415  //On the back
5416  if(s[2] < 0.0)
5417  {
5418  //It's the RUB son
5419  son = OcTreeNames::RUB;
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];
5423  }
5424  //On the front
5425  else
5426  {
5427  //It's the RUF son
5428  son = OcTreeNames::RUF;
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];
5432  }
5433  }
5434  }
5435 
5436  //Find the node in the neighbour's son
5437  Node* neighbour_son_node_pt =
5438  neigh_pt->son_pt(son)->object_pt()->
5439  get_node_at_local_coordinate(s_in_son);
5440 
5441  //If there is a node, return it
5442  if(neighbour_son_node_pt!=0)
5443  {
5444  //Return the pointer to the neighbouring node
5445  return neighbour_son_node_pt;
5446  }
5447  }
5448  }
5449  }
5450 
5451  // Keep searching, but only if there are further edge neighbours
5452  // Try next root edge neighbour
5453  i_root_edge_neighbour++;
5454 
5455  // Have we exhausted the search
5456  if (i_root_edge_neighbour>=nroot_edge_neighbour)
5457  {
5458  keep_searching=false;
5459  }
5460 
5461  } // End of while keep searching over all true edge neighbours
5462 
5463  } //End of loop over edges
5464 
5465  //Node not found, return null
5466  return 0;
5467 }
5468 
5469 //==================================================================
5470 /// Set the correct p-order of the element based on that of its
5471 /// father. Then construct an integration scheme of the correct order.
5472 /// If an adopted father is specified, information from this is
5473 /// used instead of using the father found from the tree.
5474 //==================================================================
5475 template<unsigned INITIAL_NNODE_1D>
5477 initial_setup(Tree* const &adopted_father_pt, const unsigned &initial_p_order)
5478 {
5479  //Storage for pointer to my father (in octree impersonation)
5480  OcTree* father_pt = 0;
5481 
5482  // Check if an adopted father has been specified
5483  if (adopted_father_pt!=0)
5484  {
5485  //Get pointer to my father (in octree impersonation)
5486  father_pt = dynamic_cast<OcTree*>(adopted_father_pt);
5487  }
5488  // Check if element is in a tree
5489  else if (Tree_pt!=0)
5490  {
5491  //Get pointer to my father (in octree impersonation)
5492  father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
5493  }
5494  else
5495  {
5496  //throw OomphLibError(
5497  // "Element not in a tree, and no adopted father has been specified!",
5498  // "PRefineableQElement<2,INITIAL_NNODE_1D>::initial_setup()",
5499  // OOMPH_EXCEPTION_LOCATION);
5500  }
5501 
5502  // Check if element has father
5503  if (father_pt!=0 || initial_p_order!=0)
5504  {
5505  if(father_pt!=0)
5506  {
5509  (father_pt->object_pt());
5510  if (father_el_pt!=0)
5511  {
5512  unsigned father_p_order = father_el_pt->p_order();
5513  // Set p-order to that of father
5514  P_order = father_p_order;
5515  }
5516  }
5517  else
5518  {
5519  P_order = initial_p_order;
5520  }
5521 
5522  // Now sort out the element...
5523  // (has p^3 nodes)
5524  unsigned new_n_node = P_order*P_order*P_order;
5525 
5526  // Allocate new space for Nodes (at the element level)
5527  this->set_n_node(new_n_node);
5528 
5529  // Set integration scheme
5530  delete this->integral_pt();
5531  switch(P_order)
5532  {
5533  case 2:
5534  this->set_integration_scheme(new GaussLobattoLegendre<3,2>);
5535  break;
5536  case 3:
5537  this->set_integration_scheme(new GaussLobattoLegendre<3,3>);
5538  break;
5539  case 4:
5540  this->set_integration_scheme(new GaussLobattoLegendre<3,4>);
5541  break;
5542  case 5:
5543  this->set_integration_scheme(new GaussLobattoLegendre<3,5>);
5544  break;
5545  case 6:
5546  this->set_integration_scheme(new GaussLobattoLegendre<3,6>);
5547  break;
5548  case 7:
5549  this->set_integration_scheme(new GaussLobattoLegendre<3,7>);
5550  break;
5551  default:
5552  std::ostringstream error_message;
5553  error_message <<"\nERROR: Exceeded maximum polynomial order for";
5554  error_message <<"\n integration scheme.\n";
5555  throw OomphLibError(error_message.str(),
5556  OOMPH_CURRENT_FUNCTION,
5557  OOMPH_EXCEPTION_LOCATION);
5558  }
5559  }
5560 }
5561 
5562 //==================================================================
5563 /// Check the father element for any required nodes which
5564 /// already exist
5565 //==================================================================
5566 template<unsigned INITIAL_NNODE_1D>
5568  Mesh*& mesh_pt,
5569  Vector<Node*>& new_node_pt)
5570 {
5571  using namespace OcTreeNames;
5572 
5573  //Get the number of 1d nodes
5574  unsigned n_p = nnode_1d();
5575 
5576  //Check whether static father_bound needs to be created
5577  if(Father_bound[n_p].nrow()==0) {setup_father_bounds();}
5578 
5579  //Pointer to my father (in quadtree impersonation)
5580  OcTree* father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
5581 
5582  // What type of son am I? Ask my quadtree representation...
5583  int son_type = Tree_pt->son_type();
5584 
5585  // Has somebody build me already? (If any nodes have not been built)
5586  if (!nodes_built())
5587  {
5588 #ifdef PARANOID
5589  if (father_pt==0)
5590  {
5591  std::string error_message =
5592  "Something fishy here: I have no father and yet \n";
5593  error_message +=
5594  "I have no nodes. Who has created me then?!\n";
5595 
5596  throw OomphLibError(error_message,
5597  "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
5598  OOMPH_EXCEPTION_LOCATION);
5599  }
5600 #endif
5601 
5602  // Return pointer to father element
5603  RefineableQElement<3>* father_el_pt
5604  = dynamic_cast<RefineableQElement<3>*>(father_pt->object_pt());
5605 
5606  // Timestepper should be the same for all nodes in father
5607  // element -- use it create timesteppers for new nodes
5608  TimeStepper* time_stepper_pt=father_el_pt->node_pt(0)->time_stepper_pt();
5609 
5610  // Number of history values (incl. present)
5611  unsigned ntstorage=time_stepper_pt->ntstorage();
5612 
5613  //// Pass pointer to time object:
5614  //this->time_pt()=father_el_pt->time_pt();
5615 
5616  Vector<double> s_lo(3);
5617  Vector<double> s_hi(3);
5618  Vector<double> s(3);
5619  Vector<double> x(3);
5620 
5621  // Setup vertex coordinates in father element:
5622  //--------------------------------------------
5623  switch(son_type)
5624  {
5625  case LDB:
5626  s_lo[0]=-1.0;
5627  s_hi[0]= 0.0;
5628  s_lo[1]=-1.0;
5629  s_hi[1]= 0.0;
5630  s_lo[2]=-1.0;
5631  s_hi[2]= 0.0;
5632  break;
5633 
5634  case LDF:
5635  s_lo[0]=-1.0;
5636  s_hi[0]= 0.0;
5637  s_lo[1]=-1.0;
5638  s_hi[1]= 0.0;
5639  s_lo[2]= 0.0;
5640  s_hi[2]= 1.0;
5641  break;
5642 
5643  case LUB:
5644  s_lo[0]=-1.0;
5645  s_hi[0]= 0.0;
5646  s_lo[1]= 0.0;
5647  s_hi[1]= 1.0;
5648  s_lo[2]=-1.0;
5649  s_hi[2]= 0.0;
5650  break;
5651 
5652  case LUF:
5653  s_lo[0]=-1.0;
5654  s_hi[0]= 0.0;
5655  s_lo[1]= 0.0;
5656  s_hi[1]= 1.0;
5657  s_lo[2]= 0.0;
5658  s_hi[2]= 1.0;
5659  break;
5660 
5661  case RDB:
5662  s_lo[0]= 0.0;
5663  s_hi[0]= 1.0;
5664  s_lo[1]=-1.0;
5665  s_hi[1]= 0.0;
5666  s_lo[2]=-1.0;
5667  s_hi[2]= 0.0;
5668  break;
5669 
5670  case RDF:
5671  s_lo[0]= 0.0;
5672  s_hi[0]= 1.0;
5673  s_lo[1]=-1.0;
5674  s_hi[1]= 0.0;
5675  s_lo[2]= 0.0;
5676  s_hi[2]= 1.0;
5677  break;
5678 
5679  case RUB:
5680  s_lo[0]= 0.0;
5681  s_hi[0]= 1.0;
5682  s_lo[1]= 0.0;
5683  s_hi[1]= 1.0;
5684  s_lo[2]=-1.0;
5685  s_hi[2]= 0.0;
5686  break;
5687 
5688  case RUF:
5689  s_lo[0]= 0.0;
5690  s_hi[0]= 1.0;
5691  s_lo[1]= 0.0;
5692  s_hi[1]= 1.0;
5693  s_lo[2]= 0.0;
5694  s_hi[2]= 1.0;
5695  break;
5696  }
5697 
5698  //// Pass macro element pointer on to sons and
5699  //// set coordinates in macro element
5700  //// hierher why can I see this?
5701  //if(father_el_pt->macro_elem_pt()!=0)
5702  // {
5703  // set_macro_elem_pt(father_el_pt->macro_elem_pt());
5704  // for(unsigned i=0;i<2;i++)
5705  // {
5706  // s_macro_ll(i)= father_el_pt->s_macro_ll(i)+
5707  // 0.5*(s_lo[i]+1.0)*(father_el_pt->s_macro_ur(i)-
5708  // father_el_pt->s_macro_ll(i));
5709  // s_macro_ur(i)= father_el_pt->s_macro_ll(i)+
5710  // 0.5*(s_hi[i]+1.0)*(father_el_pt->s_macro_ur(i)-
5711  // father_el_pt->s_macro_ll(i));
5712  // }
5713  // }
5714 
5715 
5716  // If the father element hasn't been generated yet, we're stuck...
5717  if(father_el_pt->node_pt(0)==0)
5718  {
5719  throw OomphLibError(
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);
5723  }
5724  else
5725  {
5726  unsigned jnod=0;
5727 
5728  Vector<double> s_fraction(3);
5729  // Loop over nodes in element
5730  for(unsigned i0=0;i0<n_p;i0++)
5731  {
5732  //Get the fractional position of the node in the direction of s[0]
5733  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
5734  // Local coordinate in father element
5735  s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
5736 
5737  for(unsigned i1=0;i1<n_p;i1++)
5738  {
5739  //Get the fractional position of the node in the direction of s[1]
5740  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
5741  // Local coordinate in father element
5742  s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
5743 
5744  for(unsigned i2=0;i2<n_p;i2++)
5745  {
5746  //Get the fractional position of the node in the direction of s[2]
5747  s_fraction[2] = this->local_one_d_fraction_of_node(i2,2);
5748  // Local coordinate in father element
5749  s[2] = s_lo[2] + (s_hi[2]-s_lo[2])*s_fraction[2];
5750 
5751  // Local node number
5752  jnod= i0 + n_p*i1 + n_p*n_p*i2;
5753 
5754  //Get the pointer to the node in the father, returns NULL
5755  //if there is not node
5756  Node* created_node_pt = father_el_pt->get_node_at_local_coordinate(s);
5757 
5758  // Does this node already exist in father element?
5759  //------------------------------------------------
5760  if(created_node_pt!=0)
5761  {
5762  // Copy node across
5763  this->node_pt(jnod) = created_node_pt;
5764 
5765  //Make sure that we update the values at the node so that
5766  //they are consistent with the present representation.
5767  //This is only need for mixed interpolation where the value
5768  //at the father could now become active.
5769 
5770  // Loop over all history values
5771  for(unsigned t=0;t<ntstorage;t++)
5772  {
5773  // Get values from father element
5774  // Note: get_interpolated_values() sets Vector size itself.
5775  Vector<double> prev_values;
5776  father_el_pt->get_interpolated_values(t,s,prev_values);
5777  //Find the minimum number of values
5778  //(either those stored at the node, or those returned by
5779  // the function)
5780  unsigned n_val_at_node = created_node_pt->nvalue();
5781  unsigned n_val_from_function = prev_values.size();
5782  //Use the ternary conditional operator here
5783  unsigned n_var = n_val_at_node < n_val_from_function ?
5784  n_val_at_node : n_val_from_function;
5785  //Assign the values that we can
5786  for(unsigned k=0;k<n_var;k++)
5787  {
5788  created_node_pt->set_value(t,k,prev_values[k]);
5789  }
5790  }
5791 
5792  }
5793  } // End of loop i2
5794  } // End of loop i1
5795  } // End of loop i0
5796  } // Sanity check: Father element has been generated
5797 
5798  } // End of nodes not built
5799 }
5800 
5801 //==================================================================
5802 /// p-refine the element inc times. (If inc<0 then p-unrefinement
5803 /// is performed.)
5804 //==================================================================
5805 template<unsigned INITIAL_NNODE_1D>
5807  const int &inc,
5808  Mesh* const &mesh_pt,
5809  GeneralisedElement* const &clone_pt)
5810 {
5811  //Cast clone to correct type
5813  = dynamic_cast<PRefineableQElement<3,INITIAL_NNODE_1D>*>(clone_pt);
5814 
5815  //Check if we can use it
5816  if(clone_el_pt==0)
5817  {
5818  throw OomphLibError(
5819  "Cloned copy must be a PRefineableQElement<3,INITIAL_NNODE_1D>!",
5820  OOMPH_CURRENT_FUNCTION,
5821  OOMPH_EXCEPTION_LOCATION);
5822  }
5823 #ifdef PARANOID
5824  //Clone exists, so check that it is infact a copy of me
5825  else
5826  {
5827  //Flag to keep track of check
5828  bool clone_is_ok = true;
5829 
5830  //Does the clone have the correct p-order?
5831  clone_is_ok = clone_is_ok && (clone_el_pt->p_order() == this->p_order());
5832 
5833  if(!clone_is_ok)
5834  {
5835  std::ostringstream err_stream;
5836  err_stream << "Clone element has a different p-order from me,"
5837  << std::endl
5838  << "but it is supposed to be a copy!" << std::endl;
5839 
5840  throw OomphLibError(err_stream.str(),
5841  OOMPH_CURRENT_FUNCTION,
5842  OOMPH_EXCEPTION_LOCATION);
5843  }
5844 
5845  //Does the clone have the same number of nodes as me?
5846  clone_is_ok = clone_is_ok && (clone_el_pt->nnode() == this->nnode());
5847 
5848  if(!clone_is_ok)
5849  {
5850  std::ostringstream err_stream;
5851  err_stream << "Clone element has a different number of nodes from me,"
5852  << std::endl
5853  << "but it is supposed to be a copy!" << std::endl;
5854 
5855  throw OomphLibError(err_stream.str(),
5856  OOMPH_CURRENT_FUNCTION,
5857  OOMPH_EXCEPTION_LOCATION);
5858  }
5859 
5860  //Does the clone have the same nodes as me?
5861  for(unsigned n=0; n<this->nnode(); n++)
5862  {
5863  clone_is_ok = clone_is_ok && (clone_el_pt->node_pt(n) == this->node_pt(n));
5864  }
5865 
5866  if(!clone_is_ok)
5867  {
5868  std::ostringstream err_stream;
5869  err_stream << "The nodes belonging to the clone element are different"
5870  << std::endl
5871  << "from mine, but it is supposed to be a copy!" << std::endl;
5872 
5873  throw OomphLibError(err_stream.str(),
5874  OOMPH_CURRENT_FUNCTION,
5875  OOMPH_EXCEPTION_LOCATION);
5876  }
5877 
5878  //If we get to here then the clone has all the information we require
5879  }
5880 #endif
5881 
5882  // Currently we can't handle the case of generalised coordinates
5883  // since we haven't established how they should be interpolated.
5884  // Buffer this case:
5885  if (clone_el_pt->node_pt(0)->nposition_type()!=1)
5886  {
5887  throw OomphLibError(
5888  "Can't handle generalised nodal positions (yet).",
5889  OOMPH_CURRENT_FUNCTION,
5890  OOMPH_EXCEPTION_LOCATION);
5891  }
5892 
5893  // Timestepper should be the same for all nodes -- use it
5894  // to create timesteppers for new nodes
5895  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
5896 
5897  // Get number of history values (incl. present)
5898  unsigned ntstorage = time_stepper_pt->ntstorage();
5899 
5900  // Increment p-order of the element
5901  this->P_order += inc;
5902 
5903  // Change integration scheme
5904  delete this->integral_pt();
5905  switch(this->P_order)
5906  {
5907  case 2:
5908  this->set_integration_scheme(new GaussLobattoLegendre<3,2>);
5909  break;
5910  case 3:
5911  this->set_integration_scheme(new GaussLobattoLegendre<3,3>);
5912  break;
5913  case 4:
5914  this->set_integration_scheme(new GaussLobattoLegendre<3,4>);
5915  break;
5916  case 5:
5917  this->set_integration_scheme(new GaussLobattoLegendre<3,5>);
5918  break;
5919  case 6:
5920  this->set_integration_scheme(new GaussLobattoLegendre<3,6>);
5921  break;
5922  case 7:
5923  this->set_integration_scheme(new GaussLobattoLegendre<3,7>);
5924  break;
5925  default:
5926  std::ostringstream error_message;
5927  error_message <<"\nERROR: Exceeded maximum polynomial order for";
5928  error_message <<"\n integration scheme.\n";
5929  throw OomphLibError(error_message.str(),
5930  OOMPH_CURRENT_FUNCTION,
5931  OOMPH_EXCEPTION_LOCATION);
5932  }
5933 
5934  // Allocate new space for Nodes (at the element level)
5935  this->set_n_node(this->P_order*this->P_order*this->P_order);
5936 
5937  // Copy vertex nodes and create new edge and internal nodes
5938  //---------------------------------------------------------
5939 
5940  // Setup vertex coordinates in element:
5941  //-------------------------------------
5942  Vector<double> s_lo(3);
5943  Vector<double> s_hi(3);
5944  s_lo[0]=-1.0;
5945  s_hi[0]= 1.0;
5946  s_lo[1]=-1.0;
5947  s_hi[1]= 1.0;
5948  s_lo[2]=-1.0;
5949  s_hi[2]= 1.0;
5950 
5951  //Local coordinate in element
5952  Vector<double> s(3);
5953 
5954  unsigned jnod=0;
5955 
5956  Vector<double> s_fraction(3);
5957  // Loop over nodes in element
5958  for(unsigned i0=0;i0<this->P_order;i0++)
5959  {
5960  //Get the fractional position of the node in the direction of s[0]
5961  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
5962  // Local coordinate
5963  s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
5964 
5965  for(unsigned i1=0;i1<this->P_order;i1++)
5966  {
5967  //Get the fractional position of the node in the direction of s[1]
5968  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
5969  // Local coordinate
5970  s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
5971 
5972  for(unsigned i2=0;i2<this->P_order;i2++)
5973  {
5974  //Get the fractional position of the node in the direction of s[2]
5975  s_fraction[2] = this->local_one_d_fraction_of_node(i2,2);
5976  // Local coordinate
5977  s[2] = s_lo[2] + (s_hi[2]-s_lo[2])*s_fraction[2];
5978 
5979  // Local node number
5980  jnod= i0 + this->P_order*i1 + this->P_order*this->P_order*i2;
5981 
5982  // Initialise flag: So far, this node hasn't been built
5983  // or copied yet
5984  bool node_done=false;
5985 
5986  //Get the pointer to the node in this element (or rather, its clone),
5987  //returns NULL if there is not node
5988  Node* created_node_pt = clone_el_pt->get_node_at_local_coordinate(s);
5989  //created_node_pt = 0;
5990 
5991  // Does this node already exist in this element?
5992  //----------------------------------------------
5993  if (created_node_pt!=0)
5994  {
5995  // Copy node across
5996  this->node_pt(jnod) = created_node_pt;
5997 
5998  //Make sure that we update the values at the node so that
5999  //they are consistent with the present representation.
6000  //This is only need for mixed interpolation where the value
6001  //at the father could now become active.
6002 
6003  // Loop over all history values
6004  for(unsigned t=0;t<ntstorage;t++)
6005  {
6006  // Get values from father element
6007  // Note: get_interpolated_values() sets Vector size itself.
6008  Vector<double> prev_values;
6009  clone_el_pt->get_interpolated_values(t,s,prev_values);
6010  //Find the minimum number of values
6011  //(either those stored at the node, or those returned by
6012  // the function)
6013  unsigned n_val_at_node = created_node_pt->nvalue();
6014  unsigned n_val_from_function = prev_values.size();
6015  //Use the ternary conditional operator here
6016  unsigned n_var = n_val_at_node < n_val_from_function ?
6017  n_val_at_node : n_val_from_function;
6018  //Assign the values that we can
6019  for(unsigned k=0;k<n_var;k++)
6020  {
6021  created_node_pt->set_value(t,k,prev_values[k]);
6022  }
6023  }
6024 
6025  // Node has been created by copy
6026  node_done = true;
6027  }
6028  // Node does not exist in this element but might already
6029  //------------------------------------------------------
6030  // have been created by neighbouring elements
6031  //-------------------------------------------
6032  else
6033  {
6034  //Was the node created by one of its neighbours
6035  //Whether or not the node lies on an edge can be calculated
6036  //by from the fractional position
6037  created_node_pt =
6038  this->node_created_by_neighbour(s_fraction);
6039 
6040  //If the node was so created, assign the pointers
6041  if (created_node_pt!=0)
6042  {
6043  this->node_pt(jnod) = created_node_pt;
6044  //Node has been created
6045  node_done = true;
6046  }
6047  // Node does not exist in neighbour element but might already
6048  //-----------------------------------------------------------
6049  // have been created by a son of a neighbouring element
6050  //-----------------------------------------------------
6051  else
6052  {
6053  //Was the node created by one of its neighbours' sons
6054  //Whether or not the node lies on an edge can be calculated
6055  //by from the fractional position
6056  created_node_pt =
6057  this->node_created_by_son_of_neighbour(s_fraction);
6058 
6059  //If the node was so created, assign the pointers
6060  if (created_node_pt!=0)
6061  {
6062  this->node_pt(jnod) = created_node_pt;
6063  //Node has been created
6064  node_done = true;
6065  } //Node does not exist in son of neighbouring element
6066  } //Node does not exist in neighbouring element
6067  } // Node does not exist in this element
6068 
6069  // Node has not been built anywhere ---> build it here
6070  if (!node_done)
6071  {
6072  //Firstly, we need to determine whether or not a node lies
6073  //on the boundary before building it, because
6074  //we actually assign a different type of node on boundaries.
6075 
6076  //Initially none
6077  int my_bound=Tree::OMEGA;
6078  //If we are on the left face
6079  if(i0==0) {my_bound = OcTreeNames::L;}
6080  //If we are on the right face
6081  else if(i0==this->P_order-1) {my_bound = OcTreeNames::R;}
6082 
6083  //If we are on the bottom face
6084  if(i1==0)
6085  {
6086  //If we already have already set the boundary, we're on an edge
6087  switch(my_bound)
6088  {
6089  case OcTreeNames::L:
6090  my_bound = OcTreeNames::LD;
6091  break;
6092  case OcTreeNames::R:
6093  my_bound = OcTreeNames::RD;
6094  break;
6095  //Boundary not set
6096  default:
6097  my_bound = OcTreeNames::D;
6098  break;
6099  }
6100  }
6101  //If we are the top face
6102  else if(i1==this->P_order-1)
6103  {
6104  //If we already have a boundary
6105  switch(my_bound)
6106  {
6107  case OcTreeNames::L:
6108  my_bound = OcTreeNames::LU;
6109  break;
6110  case OcTreeNames::R:
6111  my_bound = OcTreeNames::RU;
6112  break;
6113  default:
6114  my_bound = OcTreeNames::U;
6115  break;
6116  }
6117  }
6118 
6119  //If we are on the back face
6120  if(i2==0)
6121  {
6122  //If we already have already set the boundary, we're on an edge
6123  switch(my_bound)
6124  {
6125  case OcTreeNames::L:
6126  my_bound = OcTreeNames::LB;
6127  break;
6128  case OcTreeNames::LD:
6129  my_bound = OcTreeNames::LDB;
6130  break;
6131  case OcTreeNames::LU:
6132  my_bound = OcTreeNames::LUB;
6133  break;
6134  case OcTreeNames::R:
6135  my_bound = OcTreeNames::RB;
6136  break;
6137  case OcTreeNames::RD:
6138  my_bound = OcTreeNames::RDB;
6139  break;
6140  case OcTreeNames::RU:
6141  my_bound = OcTreeNames::RUB;
6142  break;
6143  case OcTreeNames::D:
6144  my_bound = OcTreeNames::DB;
6145  break;
6146  case OcTreeNames::U:
6147  my_bound = OcTreeNames::UB;
6148  break;
6149  //Boundary not set
6150  default:
6151  my_bound = OcTreeNames::B;
6152  break;
6153  }
6154  }
6155  //If we are the front face
6156  else if(i2==this->P_order-1)
6157  {
6158  //If we already have a boundary
6159  switch(my_bound)
6160  {
6161  case OcTreeNames::L:
6162  my_bound = OcTreeNames::LF;
6163  break;
6164  case OcTreeNames::LD:
6165  my_bound = OcTreeNames::LDF;
6166  break;
6167  case OcTreeNames::LU:
6168  my_bound = OcTreeNames::LUF;
6169  break;
6170  case OcTreeNames::R:
6171  my_bound = OcTreeNames::RF;
6172  break;
6173  case OcTreeNames::RD:
6174  my_bound = OcTreeNames::RDF;
6175  break;
6176  case OcTreeNames::RU:
6177  my_bound = OcTreeNames::RUF;
6178  break;
6179  case OcTreeNames::D:
6180  my_bound = OcTreeNames::DF;
6181  break;
6182  case OcTreeNames::U:
6183  my_bound = OcTreeNames::UF;
6184  break;
6185  default:
6186  my_bound = OcTreeNames::F;
6187  break;
6188  }
6189  }
6190 
6191  // Storage for the set of Mesh boundaries on which the
6192  // appropriate edge lives.
6193  // [New nodes should always be mid-edge nodes and therefore
6194  //only live on one boundary but just to play it safe...]
6195  std::set<unsigned> boundaries;
6196  //Only get the boundaries if we are at the edge of
6197  //an element. Nodes in the centre of an element cannot be
6198  //on Mesh boundaries
6199  if(my_bound!=Tree::OMEGA)
6200  clone_el_pt->get_boundaries(my_bound,boundaries);
6201 
6202 #ifdef PARANOID
6203  //Case where a new node lives on more than two boundaries
6204  // seems fishy enough to flag
6205  if (boundaries.size()>2)
6206  {
6207  throw OomphLibError(
6208  "boundaries.size()>2 seems a bit strange..\n",
6209  "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6210  OOMPH_EXCEPTION_LOCATION);
6211  }
6212 #endif
6213 
6214  //If the node lives on a mesh boundary,
6215  //then we need to create a boundary node
6216  if(boundaries.size()> 0)
6217  {
6218  // Create node and set the pointer to it from the element
6219  created_node_pt = this->construct_boundary_node(jnod,time_stepper_pt);
6220 
6221  //Now we need to work out whether to pin the values at
6222  //the new node based on the boundary conditions applied at
6223  //its Mesh boundary
6224 
6225  //Get the boundary conditions from the father
6226  Vector<int> bound_cons(this->ncont_interpolated_values());
6227  clone_el_pt->get_bcs(my_bound,bound_cons);
6228 
6229  //Loop over the values and pin, if necessary
6230  unsigned n_value = created_node_pt->nvalue();
6231  for(unsigned k=0;k<n_value;k++)
6232  {
6233  if (bound_cons[k]) {created_node_pt->pin(k);}
6234  }
6235 
6236  // Solid node? If so, deal with the positional boundary
6237  // conditions:
6238  SolidNode* solid_node_pt =
6239  dynamic_cast<SolidNode*>(created_node_pt);
6240  if (solid_node_pt!=0)
6241  {
6242  //Get the positional boundary conditions from the father:
6243  unsigned n_dim = created_node_pt->ndim();
6244  Vector<int> solid_bound_cons(n_dim);
6245  RefineableSolidQElement<3>* clone_solid_el_pt=
6246  dynamic_cast<RefineableSolidQElement<3>*>(clone_el_pt);
6247 #ifdef PARANOID
6248  if (clone_solid_el_pt==0)
6249  {
6250  std::string error_message =
6251  "We have a SolidNode outside a refineable SolidElement\n";
6252  error_message +=
6253  "during mesh refinement -- this doesn't make sense";
6254 
6255  throw OomphLibError(
6256  error_message,
6257  "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6258  OOMPH_EXCEPTION_LOCATION);
6259  }
6260 #endif
6261  clone_solid_el_pt->
6262  get_solid_bcs(my_bound,solid_bound_cons);
6263 
6264  //Loop over the positions and pin, if necessary
6265  for(unsigned k=0;k<n_dim;k++)
6266  {
6267  if (solid_bound_cons[k]) {solid_node_pt->pin_position(k);}
6268  }
6269  } //End of if solid_node_pt
6270 
6271  // Next, we Update the boundary lookup schemes
6272 
6273  // Check if we need to add nodes to the mesh
6274  if(mesh_pt!=0)
6275  {
6276  //Loop over the boundaries stored in the set
6277  for(std::set<unsigned>::iterator it = boundaries.begin();
6278  it != boundaries.end(); ++it)
6279  {
6280  //Add the node to the boundary
6281  mesh_pt->add_boundary_node(*it,created_node_pt);
6282 
6283  //If we have set an intrinsic coordinate on this
6284  //mesh boundary then it must also be interpolated on
6285  //the new node
6286  //Now interpolate the intrinsic boundary coordinate
6287  if(mesh_pt->boundary_coordinate_exists(*it)==true)
6288  {
6289  Vector<double> zeta(2);
6290  clone_el_pt->interpolated_zeta_on_face(*it,
6291  my_bound,
6292  s,zeta);
6293 
6294  created_node_pt->set_coordinates_on_boundary(*it,zeta);
6295  }
6296  }
6297  }
6298  }
6299  //Otherwise the node is not on a Mesh boundary and
6300  //we create a normal "bulk" node
6301  else
6302  {
6303  // Create node and set the pointer to it from the element
6304  created_node_pt = this->construct_node(jnod,time_stepper_pt);
6305  }
6306 
6307  //Now we set the position and values at the newly created node
6308 
6309  // In the first instance use macro element or FE representation
6310  // to create past and present nodal positions.
6311  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
6312  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
6313  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
6314  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
6315  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
6316  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
6317  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
6318 
6319  // Loop over # of history values
6320  for (unsigned t=0;t<ntstorage;t++)
6321  {
6322  // Get position from father element -- this uses the macro
6323  // element representation if appropriate. If the node
6324  // turns out to be a hanging node later on, then
6325  // its position gets adjusted in line with its
6326  // hanging node interpolation.
6327  Vector<double> x_prev(3);
6328  clone_el_pt->get_x(t,s,x_prev);
6329 
6330  // Set previous positions of the new node
6331  for(unsigned i=0;i<3;i++)
6332  {
6333  created_node_pt->x(t,i) = x_prev[i];
6334  }
6335  }
6336 
6337  // Loop over all history values
6338  for (unsigned t=0;t<ntstorage;t++)
6339  {
6340  // Get values from father element
6341  // Note: get_interpolated_values() sets Vector size itself.
6342  Vector<double> prev_values;
6343  clone_el_pt->get_interpolated_values(t,s,prev_values);
6344 
6345  //Initialise the values at the new node
6346  unsigned n_value = created_node_pt->nvalue();
6347  for(unsigned k=0;k<n_value;k++)
6348  {
6349  created_node_pt->set_value(t,k,prev_values[k]);
6350  }
6351  }
6352 
6353  // Add new node to mesh (if requested)
6354  if(mesh_pt!=0)
6355  {
6356  mesh_pt->add_node_pt(created_node_pt);
6357  }
6358 
6359  AlgebraicElementBase* alg_el_pt=
6360  dynamic_cast<AlgebraicElementBase*>(this);
6361 
6362  //If we do have an algebraic element
6363  if(alg_el_pt!=0)
6364  {
6365  std::string error_message =
6366  "Have not implemented p-refinement for";
6367  error_message +=
6368  "Algebraic p-refineable elements yet\n";
6369 
6370  throw
6371  OomphLibError(error_message,
6372  "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6373  OOMPH_EXCEPTION_LOCATION);
6374  }
6375 
6376  } //End of case when we build the node ourselves
6377 
6378  // Check if the element is an algebraic element
6379  AlgebraicElementBase* alg_el_pt=
6380  dynamic_cast<AlgebraicElementBase*>(this);
6381 
6382  // If the element is an algebraic element, setup
6383  // node position (past and present) from algebraic node update
6384  // function. This over-writes previous assingments that
6385  // were made based on the macro-element/FE representation.
6386  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
6387  // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
6388  // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
6389  // INFO FOR *ALL* ROOT ELEMENTS!
6390  if (alg_el_pt!=0)
6391  {
6392  // Build algebraic node update info for new node
6393  // This sets up the node update data for all node update
6394  // functions that are shared by all nodes in the father
6395  // element
6396  alg_el_pt->setup_algebraic_node_update(this->node_pt(jnod),s,
6397  clone_el_pt);
6398  }
6399 
6400  } // End of vertical loop over nodes in element (i2)
6401 
6402  } // End of horizontal loop over nodes in element (i1)
6403 
6404  } // End of horizontal loop over nodes in element (i0)
6405 
6406 
6407  // Loop over all nodes in element again, to re-set the positions
6408  // This must be done using the new element's macro-element
6409  // representation, rather than the old version which may be
6410  // of a different p-order!
6411  for(unsigned i0=0;i0<this->P_order;i0++)
6412  {
6413  //Get the fractional position of the node in the direction of s[0]
6414  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
6415  // Local coordinate
6416  s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
6417 
6418  for(unsigned i1=0;i1<this->P_order;i1++)
6419  {
6420  //Get the fractional position of the node in the direction of s[1]
6421  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
6422  // Local coordinate
6423  s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
6424 
6425  for(unsigned i2=0;i2<this->P_order;i2++)
6426  {
6427  //Get the fractional position of the node in the direction of s[2]
6428  s_fraction[2] = this->local_one_d_fraction_of_node(i2,2);
6429  // Local coordinate
6430  s[2] = s_lo[2] + (s_hi[2]-s_lo[2])*s_fraction[2];
6431 
6432  // Local node number
6433  jnod= i0 + this->P_order*i1 + this->P_order*this->P_order*i2;
6434 
6435  // Loop over # of history values
6436  for (unsigned t=0;t<ntstorage;t++)
6437  {
6438  // Get position from father element -- this uses the macro
6439  // element representation if appropriate. If the node
6440  // turns out to be a hanging node later on, then
6441  // its position gets adjusted in line with its
6442  // hanging node interpolation.
6443  Vector<double> x_prev(3);
6444  this->get_x(t,s,x_prev);
6445 
6446  // Set previous positions of the new node
6447  for(unsigned i=0;i<3;i++)
6448  {
6449  this->node_pt(jnod)->x(t,i) = x_prev[i];
6450  }
6451  }
6452  }
6453  }
6454  }
6455 
6456 
6457  // If the element is a MacroElementNodeUpdateElement, set
6458  // the update parameters for the current element's nodes --
6459  // all this needs is the vector of (pointers to the)
6460  // geometric objects that affect the MacroElement-based
6461  // node update -- this needs to be done to set the node
6462  // update info for newly created nodes
6463  MacroElementNodeUpdateElementBase* clone_m_el_pt=dynamic_cast<
6464  MacroElementNodeUpdateElementBase*>(clone_el_pt);
6465  if (clone_m_el_pt!=0)
6466  {
6467  // Get vector of geometric objects from father (construct vector
6468  // via copy operation)
6469  Vector<GeomObject*> geom_object_pt(clone_m_el_pt->geom_object_pt());
6470 
6471  // Cast current element to MacroElementNodeUpdateElement:
6472  MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
6474 
6475 #ifdef PARANOID
6476  if (m_el_pt==0)
6477  {
6478  std::string error_message =
6479  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
6480  error_message +=
6481  "Strange -- if my clone is a MacroEleme