refineable_line_element.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Non-inline and non-templated functions for RefineableQElement<1> class
31 
32 #include<algorithm>
33 
34 // oomph-lib headers
35 #include "mesh.h"
36 #include "algebraic_elements.h"
39 
40 namespace oomph
41 {
42 
43  //==========================================================================
44  /// Setup static matrix for coincidence between son nodal points and father
45  /// boundaries:
46  ///
47  /// Father_bound[nnode_1d](nnode_son,son_type) = {L/R/OMEGA}
48  ///
49  /// so that node nnode_son in element of type son_type lies on boundary/
50  /// vertex Father_bound[nnode_1d](nnode_son,son_type) in its father
51  /// element. If the node doesn't lie on a boundary the value is OMEGA.
52  //==========================================================================
54  {
55  using namespace BinaryTreeNames;
56 
57  // Find the number of nodes along a 1D edge (which is the number of nodes
58  // in the element for a 1D element!)
59  const unsigned n_node = nnode_1d();
60 
61  // Allocate space for the boundary information
62  Father_bound[n_node].resize(n_node,2);
63 
64  // Initialise: By default points are not on the boundary
65  for(unsigned n=0;n<n_node;n++)
66  {
67  for(unsigned ison=0;ison<2;ison++)
68  {
69  Father_bound[n_node](n,ison) = Tree::OMEGA;
70  }
71  }
72 
73  // Left-hand son:
74  // --------------
75 
76  // L node (0) is the L node of the parent
77  Father_bound[n_node](0,L) = L;
78 
79  // Other boundary is in the interior
80 
81  // Right-hand son:
82  // ---------------
83 
84  // R node (n_node-1) is the R node of the parent
85  Father_bound[n_node](n_node-1,R) = R;
86 
87  // Other boundary is in the interior
88  }
89 
90  //==========================================================================
91  /// If a neighbouring element has already created a node at a position
92  /// corresponding to the local fractional position within the present
93  /// element, s_fraction, return a pointer to that node. If not, return
94  /// NULL (0). If the node is periodic the flag is_periodic will be true.
95  //==========================================================================
98  bool &is_periodic)
99  {
100  using namespace BinaryTreeNames;
101 
102  // Initialise edge of current element on which node lies
103  int edge_in_current = OMEGA;
104 
105  // Determine the edge of the current element on which the node lies
106  if(s_fraction[0]==0.0) { edge_in_current = L; }
107  if(s_fraction[0]==1.0) { edge_in_current = R; }
108 
109  // If the node does not lie on an edge then there is no neighbour:
110  // return NULL
111  if(edge_in_current==OMEGA) { return 0; }
112 
113  // Allocate storage for edge in neighbouring element
114  int edge_in_neighbour;
115 
116  // Allocate storage for difference in size between current and
117  // neighbouring element
118  int diff_level;
119 
120  // Allocate storage for local coordinate of node in neighbouring tree
121  Vector<double> s_in_neighbour(1);
122 
123  // Allocate storage for flag indicating if the node is not in the same
124  // binary tree
125  bool in_neighbouring_tree;
126 
127  // Allocate storage for the pointer to the neighbouring element
128  // (using its binary tree representation)
129  BinaryTree* neighbour_pt;
130 
131  // Find pointer to neighbouring element along the edge in question and
132  // calculate the local coordinate of the node within that element
133  // s_in_neighbour
134  neighbour_pt = binary_tree_pt()->
135  gteq_edge_neighbour(edge_in_current,s_in_neighbour,edge_in_neighbour,
136  diff_level,in_neighbouring_tree);
137 
138  // If a neighbour exists...
139  if(neighbour_pt!=0)
140  {
141  // ...check whether its nodes have been created yet
142  if(neighbour_pt->object_pt()->nodes_built())
143  {
144  // If they have, find the node in question in the neighbour
145  Node* neighbour_node_pt =
146  neighbour_pt->object_pt()->get_node_at_local_coordinate(s_in_neighbour);
147 
148  // If there is no node at this position, there is a problem, since in
149  // a 1D element (whose nodes have already been built) there should
150  // ALWAYS be a node at each edge of the element.
151  if(neighbour_node_pt==0)
152  {
153  std::string error_message =
154  "Problem: an element claims to have had its nodes built, yet\n";
155  error_message +=
156  "it is missing (a least) a node at its edge.\n";
157  throw OomphLibError(error_message,
158  OOMPH_CURRENT_FUNCTION,
159  OOMPH_EXCEPTION_LOCATION);
160  }
161  // Otherwise, carry on
162  else
163  {
164  // Now work out whether it's a periodic boundary (this is only
165  // (possible if we have moved into a neighbouring tree)
166  if(in_neighbouring_tree)
167  {
168  // Return whether the neighbour is periodic
169  is_periodic =
170  binary_tree_pt()->root_pt()->is_neighbour_periodic(edge_in_current);
171  }
172  // Return the pointer to the neighbouring node
173  return neighbour_node_pt;
174  }
175  }
176  }
177  // Node not found, return null
178  return 0;
179  }
180 
181  //==========================================================================
182  /// Build the element by doing the following:
183  /// - Give it nodal positions (by establishing the pointers to its nodes).
184  /// In the process create new nodes where required (i.e. if they don't
185  /// already exist in the father element and have not already been created
186  /// whilst building new neighbouring elements). Node building involves
187  /// the following steps:
188  /// - Get the nodal position from the father element.
189  /// - Establish the time-history of the newly created nodal point
190  /// (its coordinates and previous values) consistent with the father's
191  /// history.
192  /// - Add the new node to the mesh itself.
193  /// - Doc newly created nodes in "new_nodes.dat" stored in the directory
194  /// of the DocInfo object (only if it's open!).
195  /// - NOTE: Unlike in higher-dimensional elements, in 1D it is impossible
196  /// for newly-created nodes to be on a mesh boundary, since any boundary
197  /// nodes must exist in the initial (coarse) mesh. Therefore it is not
198  /// necessary to add any nodes to the mesh's boundary node storage
199  /// schemes, and we always create normal "bulk" nodes.
200  /// - Once the element has a full complement of nodes, excute the element-
201  /// specific further_build() (empty by default -- must be overloaded for
202  /// specific elements). This deals with any build operations that are not
203  /// included in the generic process outlined above. For instance, in
204  /// Crouzeix Raviart elements we need to initialise the internal pressure
205  /// pressure values in a manner consistent with the pressure distribution
206  /// in the father element.
207  //==========================================================================
209  Vector<Node*>& new_node_pt,
210  bool& was_already_built,
211  std::ofstream &new_nodes_file)
212  {
213  using namespace BinaryTreeNames;
214 
215  // Find the number of nodes along a 1D edge (which is the number of nodes
216  // in the element for a 1D element!)
217  const unsigned n_node = nnode_1d();
218 
219  // Check whether static father_bound needs to be created
220  if(Father_bound[n_node].nrow()==0) { setup_father_bounds(); }
221 
222  // Pointer to the current element's father (in binary tree impersonation)
223  BinaryTree* father_pt
224  = dynamic_cast<BinaryTree*>(binary_tree_pt()->father_pt());
225 
226  // What type of son is the current element?
227  // Ask its binary tree representation...
228  const int son_type = Tree_pt->son_type();
229 
230  // Has somebody built the current element already?
231  // Check this by determining whether or not the first node has been built
232 
233  // If this element has not already been built...
234  if (!nodes_built())
235  {
236 #ifdef PARANOID
237  if (father_pt==0)
238  {
239  std::string error_message =
240  "Something fishy here: I have no father and yet \n";
241  error_message +=
242  "I have no nodes. Who has created me then?!\n";
243 
244  throw OomphLibError(error_message,
245  OOMPH_CURRENT_FUNCTION,
246  OOMPH_EXCEPTION_LOCATION);
247  }
248 #endif
249 
250  // Indicate status
251  was_already_built = false;
252 
253  // Return pointer to father element
254  RefineableQElement<1>* father_el_pt
255  = dynamic_cast<RefineableQElement<1>*>(father_pt->object_pt());
256 
257  // Timestepper should be the same for all nodes in father element.
258  // Use it create timesteppers for new nodes.
259  TimeStepper* time_stepper_pt=father_el_pt->node_pt(0)->time_stepper_pt();
260 
261  // Determine number of history values (including present)
262  const unsigned ntstorage = time_stepper_pt->ntstorage();
263 
264  // Currently we can't handle the case of generalised coordinates
265  // since we haven't established how they should be interpolated.
266  // Buffer this case:
267  if (father_el_pt->node_pt(0)->nposition_type()!=1)
268  {
269  throw OomphLibError(
270  "Can't handle generalised nodal positions (yet).",
271  OOMPH_CURRENT_FUNCTION,
272  OOMPH_EXCEPTION_LOCATION);
273  }
274 
275  // Set up vertex coordinates in the father element:
276  // ------------------------------------------------
277 
278  // Allocate storage for the vertex coordinates s_left and s_right of
279  // the current element as viewed by this element's father (i.e. s_left[0]
280  // stores the local coordinate within the father element at which the
281  // node on the current element's left-hand edge is located. Likewise,
282  // s_right[0] stores the local coordinate within the father element at
283  // which the node on the current element's right-hand edge is located).
284  Vector<double> s_left(1);
285  Vector<double> s_right(1);
286 
287  // In order to set up the vertex coordinates we need to know which
288  // type of son the current element is
289  switch(son_type)
290  {
291  case L:
292  s_left[0] = -1.0;
293  s_right[0] = 0.0;
294  break;
295 
296  case R:
297  s_left[0] = 0.0;
298  s_right[0] = 1.0;
299  break;
300  }
301 
302  // Pass macro element pointer on to sons and set coordinates in
303  // macro element
304  // hierher why can I see this?
305  if(father_el_pt->Macro_elem_pt!=0)
306  {
307  set_macro_elem_pt(father_el_pt->Macro_elem_pt);
308 
309  s_macro_ll(0)= father_el_pt->s_macro_ll(0)+
310  0.5*(s_left[0]+1.0)*(father_el_pt->s_macro_ur(0)-
311  father_el_pt->s_macro_ll(0));
312  s_macro_ur(0)= father_el_pt->s_macro_ll(0)+
313  0.5*(s_right[0]+1.0)*(father_el_pt->s_macro_ur(0)-
314  father_el_pt->s_macro_ll(0));
315  }
316 
317  // If the father element hasn't been generated yet, we're stuck...
318  if(father_el_pt->node_pt(0)==0)
319  {
320  throw OomphLibError(
321  "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
322  OOMPH_CURRENT_FUNCTION,
323  OOMPH_EXCEPTION_LOCATION);
324  }
325  // Otherwise, carry on
326  else
327  {
328  // Allocate storage for the location of a node in terms of the local
329  // coordinate of the father element
330  Vector<double> s_in_father(1);
331 
332  // Allocate storage for the fractional position (in the current
333  // element) of the node in the direction of s[0]
334  Vector<double> s_fraction(1);
335 
336  // Loop over all nodes in the element
337  for(unsigned n=0;n<n_node;n++)
338  {
339  // Get the fractional position (in the current element) of the node
340  // in the direction of s[0]
341  s_fraction[0] = local_one_d_fraction_of_node(n,0);
342 
343  // Evaluate the local coordinate of the node in the father element
344  s_in_father[0] = s_left[0] + (s_right[0]-s_left[0])*s_fraction[0];
345 
346  // Initialise flag: So far, this node hasn't been built or copied yet
347  bool node_done = false;
348 
349  // Get the pointer to the node in the father (returns NULL if there
350  // is no node at this position)
351  Node* created_node_pt
352  = father_el_pt->get_node_at_local_coordinate(s_in_father);
353 
354  // Does this node already exist in father element?
355  // -----------------------------------------------
356 
357  // If it does...
358  if(created_node_pt!=0)
359  {
360  // ...copy node across
361  node_pt(n) = created_node_pt;
362 
363  // Make sure that we update the values at the node so that they
364  // are consistent with the present representation. This is only
365  // needed for mixed interpolation where the value at the father
366  // could now become active.
367 
368  // Loop over all history values
369  for(unsigned t=0;t<ntstorage;t++)
370  {
371  // Get values from father element
372  // Note: get_interpolated_values() sets Vector size itself
373  Vector<double> prev_values;
374  father_el_pt->get_interpolated_values(t,s_in_father,prev_values);
375 
376  // Find the minimum number of values (either those stored at the
377  // node, or those returned by the function)
378  unsigned n_val_at_node = created_node_pt->nvalue();
379  unsigned n_val_from_function = prev_values.size();
380 
381  // Use the ternary conditional operator here
382  unsigned n_var = n_val_at_node < n_val_from_function ?
383  n_val_at_node : n_val_from_function;
384 
385  // Assign the values that we can
386  for(unsigned k=0;k<n_var;k++)
387  {
388  created_node_pt->set_value(t,k,prev_values[k]);
389  }
390  }
391 
392  // Indicate that node has been created by copy
393  node_done = true;
394  }
395 
396  // Node does not exist in father element but might already
397  // -------------------------------------------------------
398  // have been created by neighbouring elements
399  // ------------------------------------------
400 
401  else
402  {
403  // Was the node created by one of its neighbours? Whether or not
404  // the node lies on an edge can be calculated from the fractional
405  // position.
406  bool is_periodic = false;
407  created_node_pt = node_created_by_neighbour(s_fraction,is_periodic);
408 
409  // If the node was so created...
410  if(created_node_pt!=0)
411  {
412  // ...assign the pointer
413  node_pt(n) = created_node_pt;
414 
415  // Indicate that node has been created by copy
416  node_done = true;
417 
418  // In a 1D mesh there is no way that a periodic node (which must
419  // be on a boundary) can exist without being part of the initial
420  // (coarse) mesh. Therefore issue an error message if
421  // node_created_by_neighbour(...) returns `is_periodic==true'.
422 #ifdef PARANOID
423  if(is_periodic)
424  {
425  std::string error_message =
426  "node_created_by_neighbour returns a node which claims\n";
427  error_message +=
428  "to be periodic. In a 1D mesh any periodic nodes must exist\n";
429  error_message +=
430  "in the initial (coarse) mesh.";
431 
432  throw OomphLibError(error_message,
433  OOMPH_CURRENT_FUNCTION,
434  OOMPH_EXCEPTION_LOCATION);
435  }
436 #endif
437  }
438  }
439 
440  // Node does not exist in father element or neighbouring element
441  // -------------------------------------------------------------
442 
443  // If the node has not been built anywhere ---> build it here
444  if(!node_done)
445  {
446  // In a 1D mesh any node which lies on the boundary must exist in
447  // the initial (coarse) mesh, so any newly-built nodes cannot be
448  // boundary nodes. Therefore we always create a normal "bulk" node.
449 
450  // Create node and set the pointer to it from the element
451  created_node_pt = construct_node(n,time_stepper_pt);
452 
453  // Add to vector of new nodes
454  new_node_pt.push_back(created_node_pt);
455 
456  // Now we set the position and values at the newly created node.
457  // In the first instance use macro element or FE representation
458  // to create past and present nodal positions.
459  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC ELEMENTS AS NOT
460  // ALL OF THEM NECESSARILY IMPLEMENT NONTRIVIAL NODE UPDATE
461  // FUNCTIONS. CALLING THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL
462  // LEAVE THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
463  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL NOT ASSIGN SENSIBLE
464  // INITIAL POSITIONS!)
465 
466  // Loop over history values
467  for(unsigned t=0;t<ntstorage;t++)
468  {
469  // Allocate storage for the previous position of the node
470  Vector<double> x_prev(1);
471 
472  // Get position from father element -- this uses the macro
473  // element representation if appropriate.
474  father_el_pt->get_x(t,s_in_father,x_prev);
475 
476  // Set the previous position of the new node
477  created_node_pt->x(t,0) = x_prev[0];
478 
479  // Allocate storage for the previous values at the node
480  // NOTE: the size of this vector is equal to the number of values
481  // (e.g. 3 velocity components and 1 pressure, say)
482  // associated with each node and NOT the number of history values
483  // which the node stores!
484  Vector<double> prev_values;
485 
486  // Get values from father element
487  // Note: get_interpolated_values() sets Vector size itself.
488  father_el_pt->get_interpolated_values(t,s_in_father,prev_values);
489 
490  // Determine the number of values at the new node
491  const unsigned n_value = created_node_pt->nvalue();
492 
493  // Loop over all values and set the previous values
494  for(unsigned v=0;v<n_value;v++)
495  {
496  created_node_pt->set_value(t,v,prev_values[v]);
497  }
498  } // End of loop over history values
499 
500  // Add new node to mesh
501  mesh_pt->add_node_pt(created_node_pt);
502 
503  } // End of case where we build the node ourselves
504 
505  // Check if the element is an algebraic element
506  AlgebraicElementBase* alg_el_pt=
507  dynamic_cast<AlgebraicElementBase*>(this);
508 
509  // If it is, set up node position (past and present) from algebraic
510  // node position (past and present) from algebraic node update
511  // node update function. This over-writes previous assingments that
512  // were made based on the macro-element/FE representation.
513  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE NODE IS A
514  // MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN THE SAME NODAL
515  // POSITIONS BUT WE NEED TO ADD THE REMESH INFO FOR *ALL* ROOT
516  // ELEMENTS!
517  if (alg_el_pt!=0)
518  {
519  // Build algebraic node update info for new node. This sets up
520  // the node update data for all node update functions that are
521  // shared by all nodes in the father element.
522  alg_el_pt->setup_algebraic_node_update(node_pt(n),s_in_father,
523  father_el_pt);
524  }
525 
526  // If we have built the node and we are documenting our progress,
527  // write the (hopefully consistent position) to the outputfile
528  if((!node_done) && (new_nodes_file.is_open()))
529  {
530  new_nodes_file << node_pt(n)->x(0) << std::endl;
531  }
532  } // End of loop over all nodes in element
533 
534 
535  // If the element is a MacroElementNodeUpdateElement, set the update
536  // parameters for the current element's nodes -- all this needs is
537  // the vector of (pointers to the) geometric objects that affect the
538  // MacroElement-based node update. This is the same as that in the
539  // father element
540  MacroElementNodeUpdateElementBase* father_m_el_pt
541  = dynamic_cast<MacroElementNodeUpdateElementBase*>(father_el_pt);
542  if(father_m_el_pt!=0)
543  {
544  // Get Vector of geometric objects from father (construct Vector
545  // via copy operation)
546  Vector<GeomObject*> geom_object_pt(father_m_el_pt->geom_object_pt());
547 
548  // Cast current element to MacroElementNodeUpdateElement:
550  = dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
551 
552 #ifdef PARANOID
553  if (m_el_pt==0)
554  {
555  std::string error_message =
556  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
557  error_message +=
558  "Strange -- if the father is a MacroElementNodeUpdateElement\n";
559  error_message += "the son should be too....\n";
560 
561  throw OomphLibError(error_message,
562  OOMPH_CURRENT_FUNCTION,
563  OOMPH_EXCEPTION_LOCATION);
564  }
565 #endif
566  // Build update info by passing Vector of geometric objects:
567  // This sets the current element to be the update element for all
568  // of the element's nodes. This is reversed if the element is
569  // ever un-refined in the father element's rebuild_from_sons()
570  // function which overwrites this assignment to avoid nasty
571  // segmentation faults that occur when a node tries to update
572  // itself via an element that no longer exists...
573  m_el_pt->set_node_update_info(geom_object_pt);
574  }
575 
576 #ifdef OOMPH_HAS_MPI
577  // Pass on non-halo proc ID
578  Non_halo_proc_ID=
579  tree_pt()->father_pt()->object_pt()->non_halo_proc_ID();
580 #endif
581 
582  // Is the new element an ElementWithMovingNodes?
583  ElementWithMovingNodes* aux_el_pt
584  = dynamic_cast<ElementWithMovingNodes*>(this);
585 
586  // Pass down the information re the method for the evaluation
587  // of the shape derivatives
588  if(aux_el_pt!=0)
589  {
590  ElementWithMovingNodes* aux_father_el_pt
591  = dynamic_cast<ElementWithMovingNodes*>(father_el_pt);
592 
593 #ifdef PARANOID
594  if(aux_father_el_pt==0)
595  {
596  std::string error_message =
597  "Failed to cast to ElementWithMovingNodes*\n";
598  error_message +=
599  "Strange -- if the son is a ElementWithMovingNodes\n";
600  error_message += "the father should be too....\n";
601 
602  throw OomphLibError(error_message,
603  OOMPH_CURRENT_FUNCTION,
604  OOMPH_EXCEPTION_LOCATION);
605  }
606 #endif
607 
608  //If evaluating the residuals by finite differences in the father
609  //continue to do so in the child
610  if(aux_father_el_pt
611  ->are_dresidual_dnodal_coordinates_always_evaluated_by_fd())
612  {
613  aux_el_pt->
614  enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
615  }
616 
617  aux_el_pt->method_for_shape_derivs()=
618  aux_father_el_pt->method_for_shape_derivs();
619 
620  //If bypassing the evaluation of fill_in_jacobian_from_geometric_data
621  //continue to do so
622  if(aux_father_el_pt
623  ->is_fill_in_jacobian_from_geometric_data_bypassed())
624  {
626  }
627  }
628 
629 
630  // Now do further build (if any)
631  further_build();
632 
633  } // Sanity check: Father element has been generated
634 
635  } // End of if element has not already been built
636 
637  // If element has already been built, say so
638  else
639  {
640  was_already_built = true;
641  }
642  }
643 
644  //==========================================================================
645  /// Print corner nodes, use colour (default "BLACK")
646  //==========================================================================
648  std::ostream& outfile,const std::string& colour) const
649  {
650  // Allocate storage for local coordinate s
651  Vector<double> s(1);
652 
653  // Allocate storage for global coordinate of element vertex
654  Vector<double> vertex(1);
655 
656  outfile << "ZONE I=2,J=2, C=" << colour << std::endl;
657 
658  // Left-hand vertex
659  s[0]=-1.0;
660  get_x(s,vertex);
661  outfile << vertex[0] << " " << Number << std::endl;
662 
663  // Right-hand vertex
664  s[0]=1.0;
665  get_x(s,vertex);
666  outfile << vertex[0] << " " << Number << std::endl;
667 
668  outfile << "TEXT CS = GRID, X = " << vertex[0] <<
669  ", HU = GRID, H = 0.01, AN = MIDCENTER, T =\""
670  << Number << "\"" << std::endl;
671  }
672 
673  //==========================================================================
674  /// Check inter-element continuity of
675  /// - nodal positions
676  /// - (nodally) interpolated function values
677  //==========================================================================
679  {
680  using namespace BinaryTreeNames;
681 
682  // Calculate number of nodes
683  const unsigned n_node = nnode_1d();
684 
685  // Number of timesteps (including present) for which continuity is
686  // to be checked
687  const unsigned n_time=1;
688 
689  // Initialise errors
690  max_error=0.0;
691  Vector<double> max_error_x(1,0.0);
692  double max_error_val = 0.0;
693 
694  // Initialise vector of element edges in the current element
695  Vector<int> edge_in_current(2);
696  edge_in_current[0] = L;
697  edge_in_current[1] = R;
698 
699  // Loop over both edges
700  for(unsigned edge_counter=0;edge_counter<2;edge_counter++)
701  {
702  // Allocate storage for the fractional position (in the current
703  // element) of the node in the direction of s[0]
704  Vector<double> s_fraction(1);
705 
706  // Allocate storage for the location of a node in terms of the local
707  // coordinate of the current element
708  Vector<double> s_in_current(1);
709 
710  // Allocate storage for the location of a node in terms of the local
711  // coordinate of the neighbouring element
712  Vector<double> s_in_neighbour(1);
713 
714  // Allocate storage for edge in neighbouring element
715  int edge_in_neighbour;
716 
717  // Allocate storage for difference in size between current and
718  // neighbouring element
719  int diff_level;
720 
721  // Allocate storage for flag indicating if the node is not in the same
722  // binary tree
723  bool in_neighbouring_tree;
724 
725  // Calculate the local coordinate and the local coordinate as viewed
726  // from the neighbour
727 
728  // Allocate storage for the pointer to the neighbouring element
729  // (using its binary tree representation)
730  BinaryTree* neighbour_pt;
731 
732  // Find pointer to neighbouring element along the edge in question and
733  // calculate the local coordinate of the node within that element,
734  // s_in_neighbour
735  neighbour_pt = binary_tree_pt()->
736  gteq_edge_neighbour(edge_in_current[edge_counter],s_in_neighbour,
737  edge_in_neighbour,diff_level,in_neighbouring_tree);
738 
739  // If neighbour exists and has existing nodes
740  if((neighbour_pt!=0) && (neighbour_pt->object_pt()->nodes_built()))
741  {
742  // Need to exclude periodic nodes from this check.
743  // There are only periodic nodes if we are in a neighbouring tree...
744  bool is_periodic = false;
745  if(in_neighbouring_tree)
746  {
747  // Is it periodic?
748  is_periodic = tree_pt()->root_pt()
749  ->is_neighbour_periodic(edge_in_current[edge_counter]);
750  }
751 
752  // Allocate storage for pointer to the local node
753  Node* local_node_pt = 0;
754 
755  switch(edge_counter)
756  {
757  case 0:
758  // Local fraction of node
759  s_fraction[0] = 0.0;
760  // Get pointer to local node
761  local_node_pt = node_pt(0);
762  break;
763 
764  case 1:
765  // Local fraction of node
766  s_fraction[0] = 1.0;
767  // Get pointer to local node
768  local_node_pt = node_pt(n_node-1);
769  break;
770  }
771 
772  // Evaluate the local coordinate of the node in the current element
773  s_in_current[0] = -1.0 + 2.0*s_fraction[0];
774 
775  // NOTE: We have already calculated the local coordinate of the node
776  // in the neighbouring element above when calling gteq_edge_neighbour()
777 
778  // Loop over timesteps
779  for(unsigned t=0;t<n_time;t++)
780  {
781  // Allocate storage for the nodal position in the neighbouring element
782  Vector<double> x_in_neighbour(1);
783 
784  // Get the nodal position from the neighbouring element
785  neighbour_pt->object_pt()->interpolated_x(t,s_in_neighbour,
786  x_in_neighbour);
787 
788  // Check error only if the node is NOT periodic
789  if(is_periodic==false)
790  {
791  // Find the spatial error
792  double err = std::fabs(local_node_pt->x(t,0) - x_in_neighbour[0]);
793 
794  // If it's bigger than our tolerance, say so
795  if(err>1e-9)
796  {
797  oomph_info << "errx " << err << " " << t << " "
798  << local_node_pt->x(t,0)
799  << " " << x_in_neighbour[0]<< std::endl;
800 
801  oomph_info << "at " << local_node_pt->x(0) << std::endl;
802  }
803 
804  // If it's bigger than the previous max error, it is the
805  // new max error!
806  if(err>max_error_x[0]) { max_error_x[0]=err; }
807  }
808  // Allocate storage for the nodal values in the neighbouring element
809  Vector<double> values_in_neighbour;
810 
811  // Get the values from neighbouring element. NOTE: Number of values
812  // gets set by routine (because in general we don't know how many
813  // interpolated values a certain element has.
814  neighbour_pt->object_pt()
815  ->get_interpolated_values(t,s_in_neighbour,values_in_neighbour);
816 
817  // Allocate storage for the nodal values in the current element
818  Vector<double> values_in_current;
819 
820  // Get the values in current element
821  get_interpolated_values(t,s_in_current,values_in_current);
822 
823  // Now figure out how many continuously interpolated values there are
824  const unsigned num_val
825  = neighbour_pt->object_pt()->ncont_interpolated_values();
826 
827  // Check error
828  for(unsigned ival=0;ival<num_val;ival++)
829  {
830  // Find the spatial error
831  double err = std::fabs(values_in_current[ival]
832  - values_in_neighbour[ival]);
833 
834  // If it's bigger than our tolerance, say so
835  if (err>1.0e-10)
836  {
837  oomph_info << local_node_pt->x(0) << "\n# "
838  << "erru " << err << " " << ival << " "
839  << get_node_number(local_node_pt) << " "
840  << values_in_current[ival]
841  << " " << values_in_neighbour[ival] << std::endl;
842  }
843 
844  // If it's bigger than the previous max error, it is the
845  // new max error!
846  if (err>max_error_val) { max_error_val=err; }
847  }
848  } // End of loop over timesteps
849  } // End of if neighbour exists and has existing nodes
850  } // End of loop over edges
851 
852  // Update max_error if necessary
853  max_error = max_error_x[0];
854  if(max_error_val>max_error) { max_error = max_error_val; }
855 
856  // Output max_error information
857  if (max_error>1e-9)
858  {
859  oomph_info << "\n#------------------------------------ \n#Max error " ;
860  oomph_info << max_error_x[0]
861  << " " << max_error_val << std::endl;
862  oomph_info << "#------------------------------------ \n " << std::endl;
863  }
864  }
865 
866  //==========================================================================
867  /// Static matrix for coincidence between son nodal points and father
868  /// boundaries
869  //==========================================================================
870  std::map<unsigned,DenseMatrix<int> > RefineableQElement<1>::Father_bound;
871 
872 } // End of namespace
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
A policy class that serves to establish the common interfaces for elements that contain moving nodes...
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:102
void get_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
Definition: elements.h:1833
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get all continously interpolated function values in this element as a Vector. Note: Vector sets is ow...
char t
Definition: cfortran.h:572
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void setup_algebraic_node_update(Node *&node_pt, const Vector< double > &s_father, FiniteElement *father_el_pt) const
Set up node update info for (newly created) algebraic node: I.e. work out its node update information...
OomphInfo oomph_info
e
Definition: cfortran.h:575
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:966
unsigned Number
The unsigned.
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:602
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:240
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
static char t char * s
Definition: cfortran.h:572
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3804
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
Vector< GeomObject * > & geom_object_pt()
Vector of (pointers to) geometric objects involved in node update function.
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:223
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:241
Base class for elements that allow MacroElement-based node update.
int son_type() const
Return son type.
Definition: tree.h:209
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
Vector< std::string > colour
Tecplot colours.
int & method_for_shape_derivs()
Access to method (enumerated flag) for determination of shape derivs.
void enable_bypass_fill_in_jacobian_from_geometric_data()
Bypass the call to fill_in_jacobian_from_geometric_data.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
A general mesh class.
Definition: mesh.h:74
MacroElement * Macro_elem_pt
Pointer to the element's macro element (NULL by default)
Definition: elements.h:1644