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