refineable_quad_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 #include<algorithm>
31 
32 #include "mesh.h"
33 #include "algebraic_elements.h"
36 
37 namespace oomph
38 {
39 
40 
41 //==================================================================
42 /// Setup static matrix for coincidence between son nodal points and
43 /// father boundaries:
44 ///
45 /// Father_boundd[nnode_1d](nnode_son,son_type)={SW/SE/NW/NE/S/E/N/W/OMEGA}
46 ///
47 /// so that node nnode_son in element of type son_type lies
48 /// on boundary/vertex Father_boundd[nnode_1d](nnode_son,son_type) in its
49 /// father element. If the node doesn't lie on a boundary
50 /// the value is OMEGA.
51 //==================================================================
53 {
54  using namespace QuadTreeNames;
55 
56  //Find the number of nodes along a 1D edge
57  unsigned n_p = nnode_1d();
58  //Allocate space for the boundary information
59  Father_bound[n_p].resize(n_p*n_p,4);
60 
61  //Initialise: By default points are not on the boundary
62  for(unsigned n=0;n<n_p*n_p;n++)
63  {for(unsigned ison=0;ison<4;ison++) {Father_bound[n_p](n,ison)=Tree::OMEGA;}}
64 
65  // Southwest son
66  //--------------
67  // SW node (0) is the SW node of the parent
68  Father_bound[n_p](0,SW) = SW;
69  //Southern boundary is the southern boundary of the parent
70  for (unsigned n=1;n<n_p;n++) {Father_bound[n_p](n,SW) = S;}
71  //Western boundary is the western boundary of the parent
72  for(unsigned n=1;n<n_p;n++) {Father_bound[n_p](n_p*n,SW) = W;}
73  //Other boundaries are in the interior
74 
75  // Northwest son
76  //--------------
77  //NW node (n_p*(n_p-1))is the NW node of the parent
78  Father_bound[n_p](n_p*(n_p-1),NW) = NW;
79  //Northern boundary is the northern boundary of the parent
80  for(unsigned n=1;n<n_p;n++) {Father_bound[n_p](n_p*(n_p-1)+n,NW) = N;}
81  //Western boundary is the western boundary of the parent
82  for(unsigned n=0;n<(n_p-1);n++) {Father_bound[n_p](n_p*n,NW) = W;}
83  //Other boundaries are in the interior
84 
85  // Northeast son
86  //--------------
87  //NE node (n_p*n_p-1) is the NE node of the parent
88  Father_bound[n_p](n_p*n_p-1,NE) = NE;
89  //Northern boundary is the northern boundary of the parent
90  for(unsigned n=0;n<(n_p-1);n++) {Father_bound[n_p](n_p*(n_p-1)+n,NE) = N;}
91  //Eastern boundary is the eastern boundary of the parent
92  for(unsigned n=0;n<(n_p-1);n++) {Father_bound[n_p](n_p-1 + n_p*n,NE) = E;}
93  //Other boundaries are in the interior
94 
95  // Southeast son
96  //--------------
97  //SE node (n_p-1) is the SE node of the parent
98  Father_bound[n_p](n_p-1,SE) = SE;
99  //Southern boundary is the southern boundary of the parent
100  for(unsigned n=0;n<(n_p-1);n++) {Father_bound[n_p](n,SE) = S;}
101  //Eastern boundary is the eastern boundary of the parent
102  for(unsigned n=1;n<n_p;n++) {Father_bound[n_p](n_p-1 + n_p*n,SE) = E;}
103 
104 }
105 
106 //==================================================================
107 /// Determine Vector of boundary conditions along the element's boundary
108 /// (or vertex) bound (S/W/N/E/SW/SE/NW/NE).
109 ///
110 /// This function assumes that the same boundary condition is applied
111 /// along the entire length of an element's edge (of course, the
112 /// vertices combine the boundary conditions of their two adjacent edges
113 /// in the most restrictive combination. Hence, if we're at a vertex,
114 /// we apply the most restrictive boundary condition of the
115 /// two adjacent edges. If we're on an edge (in its proper interior),
116 /// we apply the least restrictive boundary condition of all nodes
117 /// along the edge.
118 ///
119 /// Usual convention:
120 /// - bound_cons[ival]=0 if value ival on this boundary is free
121 /// - bound_cons[ival]=1 if value ival on this boundary is pinned
122 //==================================================================
124  Vector<int>& bound_cons) const
125 {
126  using namespace QuadTreeNames;
127 
128  //Max. number of nodal data values in element
129  unsigned nvalue=ncont_interpolated_values();
130  //Set up temporary vectors to hold edge boundary conditions
131  Vector<int> bound_cons1(nvalue), bound_cons2(nvalue);
132 
133  //Which boundary are we on?
134  switch(bound)
135  {
136  //If on edge, just get the bcs
137  case N:
138  case S:
139  case W:
140  case E:
141  get_edge_bcs(bound,bound_cons);
142  break;
143 
144  //Most restrictive boundary at SE corner
145  case SE:
146  get_edge_bcs(S,bound_cons1);
147  get_edge_bcs(E,bound_cons2);
148 
149  for(unsigned k=0;k<nvalue;k++)
150  {bound_cons[k]=(bound_cons1[k]||bound_cons2[k]);}
151  break;
152 
153  //Most restrictive boundary at SW corner
154  case SW:
155  get_edge_bcs(S,bound_cons1);
156  get_edge_bcs(W,bound_cons2);
157 
158  for(unsigned k=0;k<nvalue;k++)
159  {bound_cons[k]=(bound_cons1[k]||bound_cons2[k]);}
160  break;
161 
162  //Most restrictive boundary at NW corner
163  case NW:
164  get_edge_bcs(N,bound_cons1);
165  get_edge_bcs(W,bound_cons2);
166 
167  for(unsigned k=0;k<nvalue;k++)
168  {bound_cons[k]=(bound_cons1[k]||bound_cons2[k]);}
169  break;
170 
171  //Most restrictive boundary at NE corner
172  case NE:
173  get_edge_bcs(N,bound_cons1);
174  get_edge_bcs(E,bound_cons2);
175 
176  for(unsigned k=0;k<nvalue;k++)
177  {bound_cons[k]=(bound_cons1[k]||bound_cons2[k]);}
178  break;
179 
180  default:
181  throw OomphLibError("Wrong boundary",
182  OOMPH_CURRENT_FUNCTION,
183  OOMPH_EXCEPTION_LOCATION);
184  }
185 }
186 
187 //==================================================================
188 /// Determine Vector of boundary conditions along the element's
189 /// edge (S/N/W/E) -- BC is the least restrictive combination
190 /// of all the nodes on this edge
191 ///
192 /// Usual convention:
193 /// - bound_cons[ival]=0 if value ival on this boundary is free
194 /// - bound_cons[ival]=1 if value ival on this boundary is pinned
195 //==================================================================
197  Vector<int>& bound_cons) const
198 {
199  using namespace QuadTreeNames;
200 
201  //Number of nodes along 1D edge
202  unsigned n_p = nnode_1d();
203  //Left- and Right-hand nodes
204  unsigned left_node, right_node;
205 
206  //Set the left (lower) and right (upper) nodes for the edge
207  switch(edge)
208  {
209  case N:
210  left_node = n_p*(n_p-1);
211  right_node = n_p*n_p - 1;
212  break;
213 
214  case S:
215  left_node = 0;
216  right_node = n_p-1;
217  break;
218 
219  case W:
220  left_node = 0;
221  right_node = n_p*(n_p-1);
222  break;
223 
224  case E:
225  left_node = n_p-1;
226  right_node = n_p*n_p - 1;
227  break;
228 
229  default:
230  std::ostringstream error_stream;
231  error_stream
232  << "Wrong edge " << edge << " passed to get_edge_bcs(..)" << std::endl;
233 
234  throw OomphLibError(error_stream.str(),
235  OOMPH_CURRENT_FUNCTION,
236  OOMPH_EXCEPTION_LOCATION);
237  }
238 
239  // Max. number of nodal data values in element
240  unsigned maxnvalue=ncont_interpolated_values();
241 
242  //Loop over all values, the least restrictive value is
243  //the boundary condition at the left multiplied by that at the right
244  //Assuming that free is always zero and pinned is one
245  for(unsigned k=0;k<maxnvalue;k++)
246  {
247  bound_cons[k] =
248  node_pt(left_node)->is_pinned(k)*node_pt(right_node)->is_pinned(k);
249  }
250 
251 }
252 
253 //==================================================================
254 /// Given an element edge/vertex, return a set that contains
255 /// all the (mesh-)boundary numbers that this element edge/vertex
256 /// lives on.
257 ///
258 /// For proper edges, the boundary is the one (if any) that is shared by
259 /// both vertex nodes). For vertex nodes, we just return their
260 /// boundaries.
261 //==================================================================
263  std::set<unsigned>& boundary) const
264 {
265  using namespace QuadTreeNames;
266 
267  //Number of 1d nodes along an edge
268  unsigned n_p = nnode_1d();
269  //Left and right-hand nodes
270  int left_node=-1, right_node=-1;
271 
272  //Set the left (lower) and right (upper) nodes for the edge
273  switch(edge)
274  {
275  case N:
276  left_node = n_p*(n_p-1);
277  right_node = n_p*n_p - 1;
278  break;
279 
280  case S:
281  left_node = 0;
282  right_node = n_p-1;
283  break;
284 
285  case W:
286  left_node = 0;
287  right_node = n_p*(n_p-1);
288  break;
289 
290  case E:
291  left_node = n_p-1;
292  right_node = n_p*n_p - 1;
293  break;
294 
295  //Vertices do not have left nodes!
296  case SE:
297  right_node = n_p-1;
298  break;
299 
300  case SW:
301  right_node = 0;
302  break;
303 
304  case NE:
305  right_node = n_p*n_p - 1;
306  break;
307 
308  case NW:
309  right_node = n_p*(n_p-1);
310  break;
311 
312  default:
313  std::ostringstream error_stream;
314  error_stream
315  << "Wrong edge " << edge << " passed" << std::endl;
316 
317  throw OomphLibError(error_stream.str(),
318  OOMPH_CURRENT_FUNCTION,
319  OOMPH_EXCEPTION_LOCATION);
320  }
321 
322  //Empty the boundary set: Edge does not live on any boundary
323  boundary.clear();
324 
325  //Storage for the nodes that the right node lives on
326  std::set<unsigned>* right_boundaries_pt=0;
327  //Get the boundaries that the right node lives on
328  node_pt(right_node)->get_boundaries_pt(right_boundaries_pt);
329 
330  //If the right node lives on some boundaries
331  if(right_boundaries_pt!=0)
332  {
333  //If the node is a vertex then add the boundaries at the node
334  //into the vector boundary
335  if(left_node < 0)
336  {
337  copy(right_boundaries_pt->begin(),
338  right_boundaries_pt->end(),
339  inserter(boundary,boundary.begin()));
340  }
341  //Otherwise only add if the boundary also exists at the left hand node
342  else
343  {
344  std::set<unsigned>* left_boundaries_pt=0;
345  node_pt(left_node)->get_boundaries_pt(left_boundaries_pt);
346  //If the left node is on some boundaries
347  if(left_boundaries_pt!=0)
348  {
349  //Use the standard algorithms to compute the boundaries in
350  //common between the left and right nodes
351  std::set_intersection(right_boundaries_pt->begin(),
352  right_boundaries_pt->end(),
353  left_boundaries_pt->begin(),
354  left_boundaries_pt->end(),
355  inserter(boundary,boundary.begin()));
356  }
357  }
358  }
359 
360 }
361 
362 
363 
364 //===================================================================
365 /// Return the value of the intrinsic boundary coordinate interpolated
366 /// along the edge (S/W/N/E)
367 //===================================================================
369 interpolated_zeta_on_edge(const unsigned &boundary,
370  const int &edge,const Vector<double> &s,
371  Vector<double> &zeta)
372 {
373  using namespace QuadTreeNames;
374 
375  //Number of 1D nodes along an edge
376  unsigned n_p = nnode_1d();
377 
378  //Storage for the shape functions
379  Shape psi(n_p*n_p);
380  //Get the shape functions at the passed position
381  this->shape(s,psi);
382 
383  //Unsigned data that give starts and multipliers for the loop
384  //over the nodes on the edges.
385  unsigned start=0, multiplier=1;
386 
387  //Which edge?
388  switch(edge)
389  {
390  case S:
391 #ifdef PARANOID
392  if(s[1] != -1.0)
393  {
394  std::ostringstream error_stream;
395  error_stream<< "Coordinate " << s[0] << " " << s[1]
396  << " is not on South edge\n";
397 
398  throw OomphLibError(error_stream.str(),
399  OOMPH_CURRENT_FUNCTION,
400  OOMPH_EXCEPTION_LOCATION);
401  }
402 #endif
403  //Start is zero and multiplier is one
404  break;
405 
406  case N:
407 #ifdef PARANOID
408  if(s[1] != 1.0)
409  {
410  std::ostringstream error_stream;
411  error_stream<< "Coordinate " << s[0] << " " << s[1]
412  << " is not on North edge\n";
413 
414  throw OomphLibError(error_stream.str(),
415  OOMPH_CURRENT_FUNCTION,
416  OOMPH_EXCEPTION_LOCATION);
417  }
418 #endif
419  //Start from the top left corner of the element, multiplier still one
420  start = n_p*(n_p-1);
421  break;
422 
423  case W:
424 #ifdef PARANOID
425  if(s[0] != -1.0)
426  {
427  std::ostringstream error_stream;
428  error_stream<< "Coordinate " << s[0] << " " << s[1]
429  << " is not on West edge\n";
430 
431  throw OomphLibError(error_stream.str(),
432  OOMPH_CURRENT_FUNCTION,
433  OOMPH_EXCEPTION_LOCATION);
434  }
435 #endif
436  //Loop over left-hand edge of element (start from zero)
437  multiplier = n_p;
438  break;
439 
440  case E:
441 #ifdef PARANOID
442  if(s[0] != 1.0)
443  {
444  std::ostringstream error_stream;
445  error_stream<< "Coordinate " << s[0] << " " << s[1]
446  << " is not on East edge\n";
447 
448  throw OomphLibError(error_stream.str(),
449  OOMPH_CURRENT_FUNCTION,
450  OOMPH_EXCEPTION_LOCATION);
451  }
452 #endif
453  //Start from the bottom right-hand corner
454  start = n_p-1;
455  //Loop over the right-hand edge of the element
456  multiplier = n_p;
457  break;
458 
459 
460  default:
461  std::ostringstream error_stream;
462  error_stream
463  << "Wrong edge " << edge << " passed" << std::endl;
464 
465  throw OomphLibError(error_stream.str(),
466  OOMPH_CURRENT_FUNCTION,
467  OOMPH_EXCEPTION_LOCATION);
468  }
469 
470  //Initialise the intrinsic coordinate
471  double inter_zeta = 0.0;
472  //Loop over the nodes on the edge
473  for(unsigned n=0;n<n_p;n++)
474  {
475  //Get the node number
476  unsigned node_number = start + multiplier*n;
477  //Now get the intrinsic coordinate
478  node_pt(node_number)->get_coordinates_on_boundary(boundary,zeta);
479  //Now multiply by the shape function
480  inter_zeta += zeta[0]*psi(node_number);
481  }
482 
483  //Set the value of the intrinsic coordinate
484  zeta[0] = inter_zeta;
485 
486 }
487 
488 
489 //===================================================================
490 /// If a neighbouring element has already created a node at
491 /// a position corresponding to the local fractional position within the
492 /// present element, s_fraction, return
493 /// a pointer to that node. If not, return NULL (0). If the node is
494 /// periodic the flag is_periodic will be true
495 //===================================================================
498  bool &is_periodic)
499 {
500  using namespace QuadTreeNames;
501 
502  //Calculate the edges on which the node lies
503  Vector<int> edges;
504  if(s_fraction[0]==0.0) {edges.push_back(W);}
505  if(s_fraction[0]==1.0) {edges.push_back(E);}
506  if(s_fraction[1]==0.0) {edges.push_back(S);}
507  if(s_fraction[1]==1.0) {edges.push_back(N);}
508 
509  //Find the number of edges
510  unsigned n_size = edges.size();
511  //If there are no edges, then there is no neighbour, return 0
512  if(n_size==0) {return 0;}
513 
514  Vector<unsigned> translate_s(2);
515  Vector<double> s_lo_neigh(2);
516  Vector<double> s_hi_neigh(2);
517  Vector<double> s(2);
518 
519  int neigh_edge, diff_level;
520  bool in_neighbouring_tree;
521 
522  //Loop over the edges
523  for(unsigned j=0;j<n_size;j++)
524  {
525  // Find pointer to neighbouring element along edge
526  QuadTree* neigh_pt;
527  neigh_pt = quadtree_pt()->
528  gteq_edge_neighbour(edges[j],translate_s,s_lo_neigh,s_hi_neigh,
529  neigh_edge,diff_level,in_neighbouring_tree);
530 
531  // Neighbour exists
532  if(neigh_pt!=0)
533  {
534  // Have its nodes been created yet?
535  //if(neigh_pt->object_pt()->node_pt(0)!=0)
536  if(neigh_pt->object_pt()->nodes_built())
537  {
538  //We now need to translate the nodal location
539  //as defined in terms of the fractional coordinates of the present
540  //element into those of its neighbour
541 
542  //Calculate the local coordinate in the neighbour
543  //Note that we need to use the translation scheme in case
544  //the local coordinates are swapped in the neighbour.
545  for(unsigned i=0;i<2;i++)
546  {
547  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
548  (s_hi_neigh[i] - s_lo_neigh[i]);
549  }
550 
551  //Find the node in the neighbour
552  Node* neighbour_node_pt =
553  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
554 
555  //If there is a node, return it
556  if(neighbour_node_pt!=0)
557  {
558  //Now work out whether it's a periodic boundary
559  //only possible if we have moved into a neighbouring tree
560  if(in_neighbouring_tree)
561  {
562  //Return whether the neighbour is periodic
563  is_periodic =
564  quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
565  }
566  //Return the pointer to the neighbouring node
567  return neighbour_node_pt;
568  }
569  }
570  }
571  }
572  //Node not found, return null
573  return 0;
574 }
575 
576 //==================================================================
577 /// Build the element by doing the following:
578 /// - Give it nodal positions (by establishing the pointers to its
579 /// nodes)
580 /// - In the process create new nodes where required (i.e. if they
581 /// don't exist in father element or have already been created
582 /// while building new neighbour elements). Node building
583 /// involves the following steps:
584 /// - Get nodal position from father element.
585 /// - Establish the time-history of the newly created nodal point
586 /// (its coordinates and the previous values) consistent with
587 /// the father's history.
588 /// - Determine the boundary conditions of the nodes (newly
589 /// created nodes can only lie on the interior of any
590 /// edges of the father element -- this makes it possible to
591 /// to figure out what their bc should be...)
592 /// - Add node to the mesh's stoarge scheme for the boundary nodes.
593 /// - Add the new node to the mesh itself
594 /// - Doc newly created nodes in "new_nodes.dat" stored in the directory
595 /// of the DocInfo object (only if it's open!)
596 /// - Finally, excute the element-specific further_build()
597 /// (empty by default -- must be overloaded for specific elements).
598 /// This deals with any build operations that are not included
599 /// in the generic process outlined above. For instance, in
600 /// Crouzeix Raviart elements we need to initialise the internal
601 /// pressure values in manner consistent with the pressure
602 /// distribution in the father element.
603 //==================================================================
605  Vector<Node*>& new_node_pt,
606  bool& was_already_built,
607  std::ofstream &new_nodes_file)
608 {
609  using namespace QuadTreeNames;
610 
611  //Get the number of 1d nodes
612  unsigned n_p = nnode_1d();
613 
614  //Check whether static father_bound needs to be created
615  if(Father_bound[n_p].nrow()==0) {setup_father_bounds();}
616 
617  //Pointer to my father (in quadtree impersonation)
618  QuadTree* father_pt = dynamic_cast<QuadTree*>(quadtree_pt()->father_pt());
619 
620  // What type of son am I? Ask my quadtree representation...
621  int son_type = Tree_pt->son_type();
622 
623  // Has somebody build me already? (If any nodes have not been built)
624  if (!nodes_built())
625  {
626 #ifdef PARANOID
627  if (father_pt==0)
628  {
629  std::string error_message =
630  "Something fishy here: I have no father and yet \n";
631  error_message +=
632  "I have no nodes. Who has created me then?!\n";
633 
634  throw OomphLibError(error_message,
635  OOMPH_CURRENT_FUNCTION,
636  OOMPH_EXCEPTION_LOCATION);
637  }
638 #endif
639 
640  // Indicate status:
641  was_already_built=false;
642 
643  // Return pointer to father element
644  RefineableQElement<2>* father_el_pt
645  = dynamic_cast<RefineableQElement<2>*>(father_pt->object_pt());
646 
647  // Timestepper should be the same for all nodes in father
648  // element -- use it create timesteppers for new nodes
649  TimeStepper* time_stepper_pt=father_el_pt->node_pt(0)->time_stepper_pt();
650 
651  // Number of history values (incl. present)
652  unsigned ntstorage=time_stepper_pt->ntstorage();
653 
654  // Currently we can't handle the case of generalised coordinates
655  // since we haven't established how they should be interpolated
656  // Buffer this case:
657  if (father_el_pt->node_pt(0)->nposition_type()!=1)
658  {
659  throw OomphLibError(
660  "Can't handle generalised nodal positions (yet).",
661  OOMPH_CURRENT_FUNCTION,
662  OOMPH_EXCEPTION_LOCATION);
663  }
664 
665  Vector<double> s_lo(2);
666  Vector<double> s_hi(2);
667  Vector<double> s(2);
668  Vector<double> x(2);
669 
670  // Setup vertex coordinates in father element:
671  //--------------------------------------------
672  switch(son_type)
673  {
674  case SW:
675  s_lo[0]=-1.0;
676  s_hi[0]= 0.0;
677  s_lo[1]=-1.0;
678  s_hi[1]= 0.0;
679  break;
680 
681  case SE:
682  s_lo[0]= 0.0;
683  s_hi[0]= 1.0;
684  s_lo[1]=-1.0;
685  s_hi[1]= 0.0;
686  break;
687 
688  case NE:
689  s_lo[0]= 0.0;
690  s_hi[0]= 1.0;
691  s_lo[1]= 0.0;
692  s_hi[1]= 1.0;
693  break;
694 
695  case NW:
696  s_lo[0]=-1.0;
697  s_hi[0]= 0.0;
698  s_lo[1]= 0.0;
699  s_hi[1]= 1.0;
700  break;
701  }
702 
703  // Pass macro element pointer on to sons and
704  // set coordinates in macro element
705  // hierher why can I see this?
706  if(father_el_pt->Macro_elem_pt!=0)
707  {
708  set_macro_elem_pt(father_el_pt->Macro_elem_pt);
709  for(unsigned i=0;i<2;i++)
710  {
711  s_macro_ll(i)= father_el_pt->s_macro_ll(i)+
712  0.5*(s_lo[i]+1.0)*(father_el_pt->s_macro_ur(i)-
713  father_el_pt->s_macro_ll(i));
714  s_macro_ur(i)= father_el_pt->s_macro_ll(i)+
715  0.5*(s_hi[i]+1.0)*(father_el_pt->s_macro_ur(i)-
716  father_el_pt->s_macro_ll(i));
717  }
718  }
719 
720 
721  // If the father element hasn't been generated yet, we're stuck...
722  if(father_el_pt->node_pt(0)==0)
723  {
724  throw OomphLibError(
725  "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
726  OOMPH_CURRENT_FUNCTION,
727  OOMPH_EXCEPTION_LOCATION);
728  }
729  else
730  {
731  unsigned jnod=0;
732  Vector<double> x_small(2);
733  Vector<double> x_large(2);
734 
735  Vector<double> s_fraction(2);
736  // Loop over nodes in element
737  for(unsigned i0=0;i0<n_p;i0++)
738  {
739  //Get the fractional position of the node in the direction of s[0]
740  s_fraction[0] = local_one_d_fraction_of_node(i0,0);
741  // Local coordinate in father element
742  s[0] = s_lo[0] + (s_hi[0]-s_lo[0])*s_fraction[0];
743 
744  for(unsigned i1=0;i1<n_p;i1++)
745  {
746  //Get the fractional position of the node in the direction of s[1]
747  s_fraction[1] = local_one_d_fraction_of_node(i1,1);
748  // Local coordinate in father element
749  s[1] = s_lo[1] + (s_hi[1]-s_lo[1])*s_fraction[1];
750 
751  // Local node number
752  jnod= i0 + n_p*i1;
753 
754  //Check whether the father's node is periodic if so, complain
755  /* {
756  Node* father_node_pt = father_el_pt->node_pt(jnod);
757  if((father_node_pt->is_a_copy()) ||
758  (father_node_pt->position_is_a_copy()))
759  {
760  throw OomphLibError(
761  "Can't handle periodic nodes (yet).",
762  OOMPH_CURRENT_FUNCTION,
763  OOMPH_EXCEPTION_LOCATION);
764  }
765  }*/
766 
767  // Initialise flag: So far, this node hasn't been built
768  // or copied yet
769  bool node_done=false;
770 
771  //Get the pointer to the node in the father, returns NULL
772  //if there is not node
773  Node* created_node_pt = father_el_pt->get_node_at_local_coordinate(s);
774 
775  // Does this node already exist in father element?
776  //------------------------------------------------
777  if(created_node_pt!=0)
778  {
779  // Copy node across
780  node_pt(jnod) = created_node_pt;
781 
782  //Make sure that we update the values at the node so that
783  //they are consistent with the present representation.
784  //This is only need for mixed interpolation where the value
785  //at the father could now become active.
786 
787  // Loop over all history values
788  for(unsigned t=0;t<ntstorage;t++)
789  {
790  // Get values from father element
791  // Note: get_interpolated_values() sets Vector size itself.
792  Vector<double> prev_values;
793  father_el_pt->get_interpolated_values(t,s,prev_values);
794  //Find the minimum number of values
795  //(either those stored at the node, or those returned by
796  // the function)
797  unsigned n_val_at_node = created_node_pt->nvalue();
798  unsigned n_val_from_function = prev_values.size();
799  //Use the ternary conditional operator here
800  unsigned n_var = n_val_at_node < n_val_from_function ?
801  n_val_at_node : n_val_from_function;
802  //Assign the values that we can
803  for(unsigned k=0;k<n_var;k++)
804  {
805  created_node_pt->set_value(t,k,prev_values[k]);
806  }
807  }
808 
809  // Node has been created by copy
810  node_done=true;
811  }
812  // Node does not exist in father element but might already
813  //--------------------------------------------------------
814  // have been created by neighbouring elements
815  //-------------------------------------------
816  else
817  {
818  //Was the node created by one of its neighbours
819  //Whether or not the node lies on an edge can be calculated
820  //by from the fractional position
821  bool is_periodic=false;;
822  created_node_pt =
823  node_created_by_neighbour(s_fraction,is_periodic);
824 
825  //If the node was so created, assign the pointers
826  if(created_node_pt!=0)
827  {
828  //If the node is periodic
829  if(is_periodic)
830  {
831  //Now the node must be on a boundary, but we don't know which
832  //one
833  //The returned created_node_pt is actually the neighbouring
834  //periodic node
835  Node* neighbour_node_pt = created_node_pt;
836 
837  // Determine the edge on which the new node will live
838  int father_bound=Father_bound[n_p](jnod,son_type);
839 
840  // Storage for the set of Mesh boundaries on which the
841  // appropriate father edge lives.
842  // [New nodes should always be mid-edge nodes in father
843  // and therefore only live on one boundary but just to
844  // play it safe...]
845  std::set<unsigned> boundaries;
846  //Only get the boundaries if we are at the edge of
847  //an element. Nodes in the centre of an element cannot be
848  //on Mesh boundaries
849  if(father_bound!=Tree::OMEGA)
850  {father_el_pt->get_boundaries(father_bound,boundaries);}
851 
852 #ifdef PARANOID
853  //Case where a new node lives on more than one boundary
854  // seems fishy enough to flag
855  if (boundaries.size()>1)
856  {
857  throw OomphLibError(
858  "boundaries.size()!=1 seems a bit strange..\n",
859  OOMPH_CURRENT_FUNCTION,
860  OOMPH_EXCEPTION_LOCATION);
861  }
862 
863  //Case when there are no boundaries, we are in big trouble
864  if(boundaries.size() == 0)
865  {
866  std::ostringstream error_stream;
867  error_stream
868  << "Periodic node is not on a boundary...\n"
869  << "Coordinates: "
870  << created_node_pt->x(0) << " "
871  << created_node_pt->x(1) << "\n";
872  throw OomphLibError(
873  error_stream.str(),
874  OOMPH_CURRENT_FUNCTION,
875  OOMPH_EXCEPTION_LOCATION);
876  }
877 #endif
878 
879  // Create node and set the pointer to it from the element
880  created_node_pt =
881  construct_boundary_node(jnod,time_stepper_pt);
882  //Make the node periodic from the neighbour
883  created_node_pt->
884  make_periodic(neighbour_node_pt);
885  // Add to vector of new nodes
886  new_node_pt.push_back(created_node_pt);
887 
888  // Loop over # of history values
889  for (unsigned t=0;t<ntstorage;t++)
890  {
891  // Get position from father element -- this uses the macro
892  // element representation if appropriate. If the node
893  // turns out to be a hanging node later on, then
894  // its position gets adjusted in line with its
895  // hanging node interpolation.
896  Vector<double> x_prev(2);
897  father_el_pt->get_x(t,s,x_prev);
898  // Set previous positions of the new node
899  for(unsigned i=0;i<2;i++)
900  {
901  created_node_pt->x(t,i) = x_prev[i];
902  }
903  }
904 
905  // Next, we Update the boundary lookup schemes
906  //Loop over the boundaries stored in the set
907  for(std::set<unsigned>::iterator it = boundaries.begin();
908  it != boundaries.end(); ++it)
909  {
910  //Add the node to the boundary
911  mesh_pt->add_boundary_node(*it,created_node_pt);
912 
913  //If we have set an intrinsic coordinate on this
914  //mesh boundary then it must also be interpolated on
915  //the new node
916  //Now interpolate the intrinsic boundary coordinate
917  if(mesh_pt->boundary_coordinate_exists(*it)==true)
918  {
919  Vector<double> zeta(1);
920  father_el_pt->interpolated_zeta_on_edge(*it,
921  father_bound,
922  s,zeta);
923 
924  created_node_pt->set_coordinates_on_boundary(*it,zeta);
925  }
926  }
927 
928  //Make sure that we add the node to the mesh
929  mesh_pt->add_node_pt(created_node_pt);
930  } //End of periodic case
931  //Otherwise the node is not periodic, so just set the
932  //pointer to the neighbours node
933  else
934  {
935  node_pt(jnod) = created_node_pt;
936  }
937  //Node has been created
938  node_done = true;
939  }
940  // Node does not exist in neighbour element but might already
941  //-----------------------------------------------------------
942  // have been created by a son of a neighbouring element
943  //-----------------------------------------------------
944  else
945  {
946  //Was the node created by one of its neighbours' sons
947  //Whether or not the node lies on an edge can be calculated
948  //by from the fractional position
949  bool is_periodic=false;;
950  created_node_pt =
951  node_created_by_son_of_neighbour(s_fraction,is_periodic);
952 
953  //If the node was so created, assign the pointers
954  if(created_node_pt!=0)
955  {
956  //If the node is periodic
957  if(is_periodic)
958  {
959  //Now the node must be on a boundary, but we don't know which
960  //one
961  //The returned created_node_pt is actually the neighbouring
962  //periodic node
963  Node* neighbour_node_pt = created_node_pt;
964 
965  // Determine the edge on which the new node will live
966  int father_bound=Father_bound[n_p](jnod,son_type);
967 
968  // Storage for the set of Mesh boundaries on which the
969  // appropriate father edge lives.
970  // [New nodes should always be mid-edge nodes in father
971  // and therefore only live on one boundary but just to
972  // play it safe...]
973  std::set<unsigned> boundaries;
974  //Only get the boundaries if we are at the edge of
975  //an element. Nodes in the centre of an element cannot be
976  //on Mesh boundaries
977  if(father_bound!=Tree::OMEGA)
978  {father_el_pt->get_boundaries(father_bound,boundaries);}
979 
980 #ifdef PARANOID
981  //Case where a new node lives on more than one boundary
982  // seems fishy enough to flag
983  if (boundaries.size()>1)
984  {
985  throw OomphLibError(
986  "boundaries.size()!=1 seems a bit strange..\n",
987  OOMPH_CURRENT_FUNCTION,
988  OOMPH_EXCEPTION_LOCATION);
989  }
990 
991  //Case when there are no boundaries, we are in big trouble
992  if(boundaries.size() == 0)
993  {
994  std::ostringstream error_stream;
995  error_stream
996  << "Periodic node is not on a boundary...\n"
997  << "Coordinates: "
998  << created_node_pt->x(0) << " "
999  << created_node_pt->x(1) << "\n";
1000  throw OomphLibError(
1001  error_stream.str(),
1002  OOMPH_CURRENT_FUNCTION,
1003  OOMPH_EXCEPTION_LOCATION);
1004  }
1005 #endif
1006 
1007  // Create node and set the pointer to it from the element
1008  created_node_pt =
1009  construct_boundary_node(jnod,time_stepper_pt);
1010  //Make the node periodic from the neighbour
1011  created_node_pt->
1012  make_periodic(neighbour_node_pt);
1013  // Add to vector of new nodes
1014  new_node_pt.push_back(created_node_pt);
1015 
1016  // Loop over # of history values
1017  for (unsigned t=0;t<ntstorage;t++)
1018  {
1019  // Get position from father element -- this uses the macro
1020  // element representation if appropriate. If the node
1021  // turns out to be a hanging node later on, then
1022  // its position gets adjusted in line with its
1023  // hanging node interpolation.
1024  Vector<double> x_prev(2);
1025  father_el_pt->get_x(t,s,x_prev);
1026  // Set previous positions of the new node
1027  for(unsigned i=0;i<2;i++)
1028  {
1029  created_node_pt->x(t,i) = x_prev[i];
1030  }
1031  }
1032 
1033  // Next, we Update the boundary lookup schemes
1034  //Loop over the boundaries stored in the set
1035  for(std::set<unsigned>::iterator it = boundaries.begin();
1036  it != boundaries.end(); ++it)
1037  {
1038  //Add the node to the boundary
1039  mesh_pt->add_boundary_node(*it,created_node_pt);
1040 
1041  //If we have set an intrinsic coordinate on this
1042  //mesh boundary then it must also be interpolated on
1043  //the new node
1044  //Now interpolate the intrinsic boundary coordinate
1045  if(mesh_pt->boundary_coordinate_exists(*it)==true)
1046  {
1047  Vector<double> zeta(1);
1048  father_el_pt->interpolated_zeta_on_edge(*it,
1049  father_bound,
1050  s,zeta);
1051 
1052  created_node_pt->set_coordinates_on_boundary(*it,zeta);
1053  }
1054  }
1055 
1056  //Make sure that we add the node to the mesh
1057  mesh_pt->add_node_pt(created_node_pt);
1058  } //End of periodic case
1059  //Otherwise the node is not periodic, so just set the
1060  //pointer to the neighbours node
1061  else
1062  {
1063  node_pt(jnod) = created_node_pt;
1064  }
1065  //Node has been created
1066  node_done = true;
1067  } // Node does not exist in son of neighbouring element
1068  } // Node does not exist in neighbouring element
1069  } // Node does not exist in father element
1070 
1071  // Node has not been built anywhere ---> build it here
1072  if (!node_done)
1073  {
1074  //Firstly, we need to determine whether or not a node lies
1075  //on the boundary before building it, because
1076  //we actually assign a different type of node on boundaries.
1077 
1078  // The node can only be on a Mesh boundary if it
1079  // lives on an edge that is shared with an edge of its
1080  // father element; i.e. it is not created inside the father element
1081  // Determine the edge on which the new node will live
1082  int father_bound=Father_bound[n_p](jnod,son_type);
1083 
1084  // Storage for the set of Mesh boundaries on which the
1085  // appropriate father edge lives.
1086  // [New nodes should always be mid-edge nodes in father
1087  // and therefore only live on one boundary but just to
1088  // play it safe...]
1089  std::set<unsigned> boundaries;
1090  //Only get the boundaries if we are at the edge of
1091  //an element. Nodes in the centre of an element cannot be
1092  //on Mesh boundaries
1093  if(father_bound!=Tree::OMEGA)
1094  {father_el_pt->get_boundaries(father_bound,boundaries);}
1095 
1096 #ifdef PARANOID
1097  //Case where a new node lives on more than one boundary
1098  // seems fishy enough to flag
1099  if (boundaries.size()>1)
1100  {
1101  throw OomphLibError(
1102  "boundaries.size()!=1 seems a bit strange..\n",
1103  OOMPH_CURRENT_FUNCTION,
1104  OOMPH_EXCEPTION_LOCATION);
1105  }
1106 #endif
1107 
1108  //If the node lives on a mesh boundary,
1109  //then we need to create a boundary node
1110  if(boundaries.size()> 0)
1111  {
1112  // Create node and set the pointer to it from the element
1113  created_node_pt = construct_boundary_node(jnod,time_stepper_pt);
1114  // Add to vector of new nodes
1115  new_node_pt.push_back(created_node_pt);
1116 
1117  //Now we need to work out whether to pin the values at
1118  //the new node based on the boundary conditions applied at
1119  //its Mesh boundary
1120 
1121  //Get the boundary conditions from the father
1122  Vector<int> bound_cons(ncont_interpolated_values());
1123  father_el_pt->get_bcs(father_bound,bound_cons);
1124 
1125  //Loop over the values and pin, if necessary
1126  unsigned n_value = created_node_pt->nvalue();
1127  for(unsigned k=0;k<n_value;k++)
1128  {
1129  if (bound_cons[k]) {created_node_pt->pin(k);}
1130  }
1131 
1132  // Solid node? If so, deal with the positional boundary
1133  // conditions:
1134  SolidNode* solid_node_pt =
1135  dynamic_cast<SolidNode*>(created_node_pt);
1136  if (solid_node_pt!=0)
1137  {
1138  //Get the positional boundary conditions from the father:
1139  unsigned n_dim = created_node_pt->ndim();
1140  Vector<int> solid_bound_cons(n_dim);
1141  RefineableSolidQElement<2>* father_solid_el_pt=
1142  dynamic_cast<RefineableSolidQElement<2>*>(father_el_pt);
1143 #ifdef PARANOID
1144  if (father_solid_el_pt==0)
1145  {
1146  std::string error_message =
1147  "We have a SolidNode outside a refineable SolidElement\n";
1148  error_message +=
1149  "during mesh refinement -- this doesn't make sense";
1150 
1151  throw OomphLibError(error_message,
1152  OOMPH_CURRENT_FUNCTION,
1153  OOMPH_EXCEPTION_LOCATION);
1154  }
1155 #endif
1156  father_solid_el_pt->
1157  get_solid_bcs(father_bound,solid_bound_cons);
1158 
1159  //Loop over the positions and pin, if necessary
1160  for(unsigned k=0;k<n_dim;k++)
1161  {
1162  if (solid_bound_cons[k]) {solid_node_pt->pin_position(k);}
1163  }
1164  } //End of if solid_node_pt
1165 
1166 
1167 
1168  // Next, we Update the boundary lookup schemes
1169  //Loop over the boundaries stored in the set
1170  for(std::set<unsigned>::iterator it = boundaries.begin();
1171  it != boundaries.end(); ++it)
1172  {
1173  //Add the node to the boundary
1174  mesh_pt->add_boundary_node(*it,created_node_pt);
1175 
1176  //If we have set an intrinsic coordinate on this
1177  //mesh boundary then it must also be interpolated on
1178  //the new node
1179  //Now interpolate the intrinsic boundary coordinate
1180  if(mesh_pt->boundary_coordinate_exists(*it)==true)
1181  {
1182  Vector<double> zeta(1);
1183  father_el_pt->interpolated_zeta_on_edge(*it,
1184  father_bound,
1185  s,zeta);
1186 
1187  created_node_pt->set_coordinates_on_boundary(*it,zeta);
1188  }
1189  }
1190  }
1191  //Otherwise the node is not on a Mesh boundary and
1192  //we create a normal "bulk" node
1193  else
1194  {
1195  // Create node and set the pointer to it from the element
1196  created_node_pt = construct_node(jnod,time_stepper_pt);
1197  // Add to vector of new nodes
1198  new_node_pt.push_back(created_node_pt);
1199  }
1200 
1201  //Now we set the position and values at the newly created node
1202 
1203  // In the first instance use macro element or FE representation
1204  // to create past and present nodal positions.
1205  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
1206  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
1207  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
1208  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
1209  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
1210  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
1211  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
1212 
1213  // Loop over # of history values
1214  for (unsigned t=0;t<ntstorage;t++)
1215  {
1216  // Get position from father element -- this uses the macro
1217  // element representation if appropriate. If the node
1218  // turns out to be a hanging node later on, then
1219  // its position gets adjusted in line with its
1220  // hanging node interpolation.
1221  Vector<double> x_prev(2);
1222  father_el_pt->get_x(t,s,x_prev);
1223 
1224  // Set previous positions of the new node
1225  for(unsigned i=0;i<2;i++)
1226  {
1227  created_node_pt->x(t,i) = x_prev[i];
1228  }
1229  }
1230 
1231  // Loop over all history values
1232  for (unsigned t=0;t<ntstorage;t++)
1233  {
1234  // Get values from father element
1235  // Note: get_interpolated_values() sets Vector size itself.
1236  Vector<double> prev_values;
1237  father_el_pt->get_interpolated_values(t,s,prev_values);
1238  //Initialise the values at the new node
1239  unsigned n_value = created_node_pt->nvalue();
1240  for(unsigned k=0;k<n_value;k++)
1241  {
1242  created_node_pt->set_value(t,k,prev_values[k]);
1243  }
1244  }
1245 
1246  // Add new node to mesh
1247  mesh_pt->add_node_pt(created_node_pt);
1248 
1249  } //End of case when we build the node ourselves
1250 
1251  // Check if the element is an algebraic element
1252  AlgebraicElementBase* alg_el_pt=
1253  dynamic_cast<AlgebraicElementBase*>(this);
1254 
1255  // If the element is an algebraic element, setup
1256  // node position (past and present) from algebraic node update
1257  // function. This over-writes previous assingments that
1258  // were made based on the macro-element/FE representation.
1259  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
1260  // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
1261  // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
1262  // INFO FOR *ALL* ROOT ELEMENTS!
1263  if (alg_el_pt!=0)
1264  {
1265  // Build algebraic node update info for new node
1266  // This sets up the node update data for all node update
1267  // functions that are shared by all nodes in the father
1268  // element
1269  alg_el_pt->setup_algebraic_node_update(node_pt(jnod),s,
1270  father_el_pt);
1271  }
1272 
1273  // If we have built the node and we are documenting our progress
1274  // write the (hopefully consistent position) to the outputfile
1275  if((!node_done) && (new_nodes_file.is_open()))
1276  {
1277  new_nodes_file << node_pt(jnod)->x(0) << " "
1278  << node_pt(jnod)->x(1) << std::endl;
1279  }
1280 
1281  } // End of vertical loop over nodes in element
1282 
1283  } // End of horizontal loop over nodes in element
1284 
1285 
1286  // If the element is a MacroElementNodeUpdateElement, set
1287  // the update parameters for the current element's nodes --
1288  // all this needs is the vector of (pointers to the)
1289  // geometric objects that affect the MacroElement-based
1290  // node update -- this is the same as that in the father element
1291  MacroElementNodeUpdateElementBase* father_m_el_pt=dynamic_cast<
1292  MacroElementNodeUpdateElementBase*>(father_el_pt);
1293  if (father_m_el_pt!=0)
1294  {
1295  // Get vector of geometric objects from father (construct vector
1296  // via copy operation)
1297  Vector<GeomObject*> geom_object_pt(father_m_el_pt->geom_object_pt());
1298 
1299  // Cast current element to MacroElementNodeUpdateElement:
1300  MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
1302 
1303 #ifdef PARANOID
1304  if (m_el_pt==0)
1305  {
1306  std::string error_message =
1307  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
1308  error_message +=
1309  "Strange -- if the father is a MacroElementNodeUpdateElement\n";
1310  error_message += "the son should be too....\n";
1311 
1312  throw OomphLibError(error_message,
1313  OOMPH_CURRENT_FUNCTION,
1314  OOMPH_EXCEPTION_LOCATION);
1315  }
1316 #endif
1317  // Build update info by passing vector of geometric objects:
1318  // This sets the current element to be the update element
1319  // for all of the element's nodes -- this is reversed
1320  // if the element is ever un-refined in the father element's
1321  // rebuild_from_sons() function which overwrites this
1322  // assignment to avoid nasty segmentation faults that occur
1323  // when a node tries to update itself via an element that no
1324  // longer exists...
1325  m_el_pt->set_node_update_info(geom_object_pt);
1326  }
1327 
1328 #ifdef OOMPH_HAS_MPI
1329  // Pass on non-halo proc id
1330  Non_halo_proc_ID=
1331  tree_pt()->father_pt()->object_pt()->non_halo_proc_ID();
1332 #endif
1333 
1334  // Is it an ElementWithMovingNodes?
1335  ElementWithMovingNodes* aux_el_pt=
1336  dynamic_cast<ElementWithMovingNodes*>(this);
1337 
1338  // Pass down the information re the method for the evaluation
1339  // of the shape derivatives
1340  if (aux_el_pt!=0)
1341  {
1342  ElementWithMovingNodes* aux_father_el_pt=
1343  dynamic_cast<ElementWithMovingNodes*>(father_el_pt);
1344 
1345 #ifdef PARANOID
1346  if (aux_father_el_pt==0)
1347  {
1348  std::string error_message =
1349  "Failed to cast to ElementWithMovingNodes*\n";
1350  error_message +=
1351  "Strange -- if the son is a ElementWithMovingNodes\n";
1352  error_message += "the father should be too....\n";
1353 
1354  throw OomphLibError(error_message,
1355  OOMPH_CURRENT_FUNCTION,
1356  OOMPH_EXCEPTION_LOCATION);
1357  }
1358 #endif
1359 
1360  //If evaluating the residuals by finite differences in the father
1361  //continue to do so in the child
1362  if(aux_father_el_pt
1363  ->are_dresidual_dnodal_coordinates_always_evaluated_by_fd())
1364  {
1365  aux_el_pt->
1366  enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
1367  }
1368 
1369  aux_el_pt->method_for_shape_derivs()=
1370  aux_father_el_pt->method_for_shape_derivs();
1371 
1372  //If bypassing the evaluation of fill_in_jacobian_from_geometric_data
1373  //continue to do so
1374  if(aux_father_el_pt
1375  ->is_fill_in_jacobian_from_geometric_data_bypassed())
1376  {
1378  }
1379  }
1380 
1381 
1382  // Now do further build (if any)
1383  further_build();
1384 
1385  } // Sanity check: Father element has been generated
1386 
1387  }
1388  // Element has already been built
1389  else
1390  {
1391  was_already_built=true;
1392  }
1393 
1394 }
1395 
1396 //====================================================================
1397 /// Print corner nodes, use colour (default "BLACK")
1398 //====================================================================
1399 void RefineableQElement<2>::output_corners(std::ostream& outfile,
1400  const std::string& colour) const
1401 {
1402  Vector<double> s(2);
1403  Vector<double> corner(2);
1404 
1405  outfile << "ZONE I=2,J=2, C=" << colour << std::endl;
1406 
1407  s[0]=-1.0;
1408  s[1]=-1.0;
1409  get_x(s,corner);
1410  outfile << corner[0] << " " << corner[1] << " "
1411  << Number << std::endl;
1412 
1413  s[0]=1.0;
1414  s[1]=-1.0;
1415  get_x(s,corner);
1416  outfile << corner[0] << " " << corner[1] << " "
1417  << Number << std::endl;
1418 
1419  s[0]=-1.0;
1420  s[1]=1.0;
1421  get_x(s,corner);
1422  outfile << corner[0] << " " << corner[1] << " "
1423  << Number << std::endl;
1424 
1425  s[0]=1.0;
1426  s[1]=1.0;
1427  get_x(s,corner);
1428  outfile << corner[0] << " " << corner[1] << " "
1429  << Number << std::endl;
1430 
1431 
1432  outfile << "TEXT CS = GRID, X = " << corner[0] <<
1433  ",Y = " << corner[1] << ", HU = GRID, H = 0.01, AN = MIDCENTER, T =\""
1434  << Number << "\"" << std::endl;
1435 
1436 }
1437 
1438 //====================================================================
1439 /// Set up all hanging nodes. If we are documenting the output then
1440 /// open the output files and pass the open files to the helper function
1441 //====================================================================
1443  &output_stream)
1444 {
1445 #ifdef PARANOID
1446  if(output_stream.size() != 4)
1447  {
1448  throw OomphLibError("There must be four output streams",
1449  OOMPH_CURRENT_FUNCTION,
1450  OOMPH_EXCEPTION_LOCATION);
1451  }
1452 #endif
1453 
1454  using namespace QuadTreeNames;
1455 
1456  //Setup hanging nodes on each edge of the element
1457  quad_hang_helper(-1,S,*(output_stream[0]));
1458  quad_hang_helper(-1,N,*(output_stream[1]));
1459  quad_hang_helper(-1,W,*(output_stream[2]));
1460  quad_hang_helper(-1,E,*(output_stream[3]));
1461 }
1462 
1463 //================================================================
1464 /// Internal function that sets up the hanging node scheme for
1465 /// a particular continuously interpolated value
1466 //===============================================================
1468 {
1469  using namespace QuadTreeNames;
1470 
1471  std::ofstream dummy_hangfile;
1472  quad_hang_helper(value_id,S,dummy_hangfile);
1473  quad_hang_helper(value_id,N,dummy_hangfile);
1474  quad_hang_helper(value_id,W,dummy_hangfile);
1475  quad_hang_helper(value_id,E,dummy_hangfile);
1476 }
1477 
1478 
1479 //=================================================================
1480 /// Internal function to set up the hanging nodes on a particular
1481 /// edge of the element
1482 //=================================================================
1484 quad_hang_helper(const int &value_id,
1485  const int &my_edge, std::ofstream& output_hangfile)
1486 {
1487  using namespace QuadTreeNames;
1488 
1489  Vector<unsigned> translate_s(2);
1490  Vector<double> s_lo_neigh(2);
1491  Vector<double> s_hi_neigh(2);
1492  int neigh_edge,diff_level;
1493  bool in_neighbouring_tree;
1494 
1495  // Find pointer to neighbour in this direction
1496  QuadTree* neigh_pt;
1497  neigh_pt=quadtree_pt()->
1498  gteq_edge_neighbour(my_edge, translate_s, s_lo_neigh,
1499  s_hi_neigh,neigh_edge,diff_level,in_neighbouring_tree);
1500 
1501  // Neighbour exists and all nodes have been created
1502  if(neigh_pt!=0)
1503  {
1504  // Different sized element?
1505  if(diff_level!=0)
1506  {
1507  //Test for the periodic node case
1508  //Are we crossing a periodic boundary
1509  bool is_periodic = false;
1510  if(in_neighbouring_tree)
1511  {is_periodic = tree_pt()->root_pt()->is_neighbour_periodic(my_edge);}
1512 
1513  //If it is periodic we actually need to get the node in
1514  //the neighbour of the neighbour (which will be a parent of
1515  //the present element) so that the "fixed" coordinate is
1516  //correctly calculated.
1517  //The idea is to replace the neigh_pt and associated data
1518  //with those of the neighbour of the neighbour
1519  if(is_periodic)
1520  {
1521  //Required data for the neighbour finding routine
1522  Vector<unsigned> translate_s_in_neigh(2);
1523  Vector<double> s_lo_neigh_of_neigh(2);
1524  Vector<double> s_hi_neigh_of_neigh(2);
1525  int neigh_edge_of_neigh,diff_level_of_neigh;
1526  bool in_neighbouring_tree_of_neigh;
1527 
1528  // Find pointer to neighbour of the neighbour on the edge
1529  // that we are currently considering
1530  QuadTree* neigh_of_neigh_pt;
1531  neigh_of_neigh_pt = neigh_pt->
1532  gteq_edge_neighbour(neigh_edge, translate_s_in_neigh,
1533  s_lo_neigh_of_neigh,
1534  s_hi_neigh_of_neigh,
1535  neigh_edge_of_neigh,diff_level_of_neigh,
1536  in_neighbouring_tree_of_neigh);
1537 
1538  //Set the value of the NEW neighbour and edge
1539  neigh_pt = neigh_of_neigh_pt;
1540  neigh_edge = neigh_edge_of_neigh;
1541 
1542  //Set the values of the translation schemes
1543  //Need to find the values of s_lo and s_hi
1544  //in the neighbour of the neighbour
1545 
1546  //Get the minimum and maximum values of the coordinate
1547  //in the neighbour (don't like this, but I think it's
1548  //necessary) Note that these values are hardcoded
1549  //in the quadtrees at some point!!
1550  double s_min = neigh_pt->object_pt()->s_min();
1551  double s_max = neigh_pt->object_pt()->s_max();
1552  Vector<double> s_lo_frac(2), s_hi_frac(2);
1553  //Work out the fractional position of the low and high points
1554  //of the original element
1555  for(unsigned i=0;i<2;i++)
1556  {
1557  s_lo_frac[i] = (s_lo_neigh[i] - s_min)/(s_max - s_min);
1558  s_hi_frac[i] = (s_hi_neigh[i] - s_min)/(s_max - s_min);
1559  }
1560 
1561  //We should now be able to construct the low and high points in
1562  //the neighbour of the neighbour
1563  for(unsigned i=0;i<2;i++)
1564  {
1565  s_lo_neigh[i] = s_lo_neigh_of_neigh[i]
1566  + s_lo_frac[translate_s_in_neigh[i]]*
1567  (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
1568  s_hi_neigh[i] = s_lo_neigh_of_neigh[i]
1569  + s_hi_frac[translate_s_in_neigh[i]]*
1570  (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
1571  }
1572 
1573  //Finally we must sort out the translation scheme
1574  Vector<unsigned> temp_translate(2);
1575  for(unsigned i=0;i<2;i++)
1576  {
1577  temp_translate[i] = translate_s_in_neigh[translate_s[i]];
1578  }
1579  for(unsigned i=0;i<2;i++) {translate_s[i] = temp_translate[i];}
1580  } //End of special treatment for periodic hanging nodes
1581 
1582  //Number of nodes in one dimension
1583  unsigned n_p = ninterpolating_node_1d(value_id);
1584  //Storage for the local nodes along the edge of the quadtree
1585  Node* local_node_pt=0;
1586  // Loop over nodes along the edge
1587  for(unsigned i0=0;i0<n_p;i0++)
1588  {
1589  //Storage for the fractional position of the node in the element
1590  Vector<double> s_fraction(2);
1591 
1592  // Find the local node and the fractional position of the node
1593  // which depends on the edge, of course
1594  switch(my_edge)
1595  {
1596  case N:
1597  s_fraction[0] =
1598  local_one_d_fraction_of_interpolating_node(i0,0,value_id);
1599  s_fraction[1] = 1.0;
1600  local_node_pt = interpolating_node_pt(i0 + n_p*(n_p-1),value_id);
1601  break;
1602 
1603  case S:
1604  s_fraction[0] =
1605  local_one_d_fraction_of_interpolating_node(i0,0,value_id);
1606  s_fraction[1] = 0.0;
1607  local_node_pt = interpolating_node_pt(i0,value_id);
1608  break;
1609 
1610  case E:
1611  s_fraction[0] = 1.0;
1612  s_fraction[1] =
1613  local_one_d_fraction_of_interpolating_node(i0,1,value_id);
1614  local_node_pt = interpolating_node_pt(n_p-1 + n_p*i0,value_id);
1615  break;
1616 
1617  case W:
1618  s_fraction[0] = 0.0;
1619  s_fraction[1] =
1620  local_one_d_fraction_of_interpolating_node(i0,1,value_id);
1621  local_node_pt = interpolating_node_pt(n_p*i0,value_id);
1622  break;
1623 
1624  default:
1625  throw OomphLibError("my_edge not N, S, W, E\n",
1626  OOMPH_CURRENT_FUNCTION,
1627  OOMPH_EXCEPTION_LOCATION);
1628  }
1629 
1630  //Calculate the local coordinates of the node in the neighbour
1631  Vector<double> s_in_neighb(2);
1632  for(unsigned i=0;i<2;i++)
1633  {
1634  s_in_neighb[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
1635  (s_hi_neigh[i] - s_lo_neigh[i]);
1636  }
1637 
1638  //Find the Node in the neighbouring element
1639  Node* neighbouring_node_pt = neigh_pt->object_pt()
1640  ->get_interpolating_node_at_local_coordinate(s_in_neighb,value_id);
1641 
1642  //If the neighbour does not have a node at this point
1643  if(0==neighbouring_node_pt)
1644  {
1645  //Do we need to make a hanging node, we assume that we don't
1646  //initially
1647  bool make_hanging_node = false;
1648 
1649  // If the node is not hanging geometrically, then we must make
1650  // it hang
1651  if(!local_node_pt->is_hanging())
1652  {
1653  make_hanging_node = true;
1654  }
1655  //Otherwise, it could be hanging geometrically, but still
1656  //require a different hanging scheme for this data value
1657  else
1658  {
1659  if(local_node_pt->hanging_pt(value_id)==local_node_pt->hanging_pt())
1660  {
1661  make_hanging_node = true;
1662  }
1663  }
1664 
1665  //If we do need to make the hanging node, then let's do it
1666  if(make_hanging_node == true)
1667  {
1668  //Cache refineable element used here
1669  RefineableElement* const obj_pt = neigh_pt->object_pt();
1670 
1671  //Get shape functions in neighbour element
1672  Shape psi(obj_pt->ninterpolating_node(value_id));
1673  obj_pt->interpolating_basis(s_in_neighb,psi,value_id);
1674 
1675  //Allocate the storage for the Hang pointer
1676  //which contains n_p nodes
1677  HangInfo *hang_pt = new HangInfo(n_p);
1678 
1679  // Loop over nodes on edge in neighbour and mark them as nodes
1680  // that this node depends on
1681  unsigned n_neighbour;
1682 
1683  // Number of nodes along edge in neighbour element
1684  for(unsigned n_edge=0;n_edge<n_p;n_edge++)
1685  {
1686  switch(neigh_edge)
1687  {
1688  case N:
1689  n_neighbour = n_p*(n_p-1) + n_edge;
1690  break;
1691 
1692  case S:
1693  n_neighbour = n_edge;
1694  break;
1695 
1696  case W:
1697  n_neighbour = n_p*n_edge;
1698  break;
1699 
1700  case E:
1701  n_neighbour = n_p*n_edge + (n_p-1);
1702  break;
1703 
1704  default:
1705  throw OomphLibError("neigh_edge not N, S, W, E\n",
1706  OOMPH_CURRENT_FUNCTION,
1707  OOMPH_EXCEPTION_LOCATION);
1708  }
1709 
1710  // Push back neighbouring node and weight into
1711  // Vector of (pointers to)
1712  // master nodes and weights
1713  // The weight is merely the value of the shape function
1714  //corresponding to the node in the neighbour
1715  hang_pt->set_master_node_pt(
1716  n_edge,obj_pt->interpolating_node_pt(n_neighbour,value_id),
1717  psi[n_neighbour]);
1718  }
1719 
1720  //Now set the hanging data for the position
1721  //This also constrains the data values associated with the
1722  //value id
1723  local_node_pt->set_hanging_pt(hang_pt,value_id);
1724  }
1725 
1726  //Dump the output if the file has been openeed
1727  if(output_hangfile.is_open())
1728  {
1729  output_hangfile << local_node_pt->x(0) << " " <<
1730  local_node_pt->x(1) << std::endl;
1731  }
1732  }
1733  //Otherwise check that the nodes are the same
1734  else
1735  {
1736 #ifdef PARANOID
1737  if (local_node_pt!=neighbouring_node_pt)
1738  {
1739  std::ostringstream warning_stream;
1740  warning_stream << "SANITY CHECK in quad_hang_helper \n"
1741  << "Current node " << local_node_pt
1742  << " at "
1743  << "(" << local_node_pt->x(0) << ", "
1744  << local_node_pt->x(1) << ")"
1745  << " is not hanging and has " << std::endl
1746  << "Neighbour's node " << neighbouring_node_pt
1747  << " at "
1748  << "(" << neighbouring_node_pt->x(0) << ", "
1749  << neighbouring_node_pt->x(1) << ")"
1750  << std::endl
1751  << "even though the two should be "
1752  << "identical" << std::endl;
1753  OomphLibWarning(warning_stream.str(),
1754  "RefineableQElement<2>::quad_hang_helper()",
1755  OOMPH_EXCEPTION_LOCATION);
1756  }
1757 #endif
1758  }
1759 
1760  //If we are doing the position, then
1761  if(value_id==-1)
1762  {
1763  // Get the nodal position from neighbour element
1764  Vector<double> x_in_neighb(2);
1765  neigh_pt->object_pt()->interpolated_x(s_in_neighb,x_in_neighb);
1766 
1767  // Fine adjust the coordinates (macro map will pick up boundary
1768  // accurately but will lead to different element edges)
1769  local_node_pt->x(0)=x_in_neighb[0];
1770  local_node_pt->x(1)=x_in_neighb[1];
1771  }
1772  }
1773  }
1774  }
1775 }
1776 
1777 
1778 //=================================================================
1779 /// Check inter-element continuity of
1780 /// - nodal positions
1781 /// - (nodally) interpolated function values
1782 //====================================================================
1783 //template<unsigned NNODE_1D>
1785 {
1786 
1787  using namespace QuadTreeNames;
1788 
1789  // Number of nodes along edge
1790  unsigned n_p=nnode_1d();
1791 
1792  // Number of timesteps (incl. present) for which continuity is
1793  // to be checked.
1794  unsigned n_time=1;
1795 
1796  // Initialise errors
1797  max_error=0.0;
1798  Vector<double> max_error_x(2,0.0);
1799  double max_error_val=0.0;
1800 
1801  Vector<int> edges(4);
1802  edges[0] = S; edges[1] = N; edges[2] = W; edges[3] = E;
1803 
1804  //Loop over the edges
1805  for(unsigned edge_counter=0;edge_counter<4;edge_counter++)
1806  {
1807  Vector<unsigned> translate_s(2);
1808  Vector<double> s(2), s_lo_neigh(2), s_hi_neigh(2), s_fraction(2);
1809  int neigh_edge,diff_level;
1810  bool in_neighbouring_tree;
1811 
1812  // Find pointer to neighbour in this direction
1813  QuadTree* neigh_pt;
1814  neigh_pt=quadtree_pt()->gteq_edge_neighbour(edges[edge_counter],
1815  translate_s,
1816  s_lo_neigh,s_hi_neigh,
1817  neigh_edge,diff_level,
1818  in_neighbouring_tree);
1819 
1820  // Neighbour exists and has existing nodes
1821  if((neigh_pt!=0) && (neigh_pt->object_pt()->nodes_built()))
1822  {
1823  //Need to exclude periodic nodes from this check
1824  //There are only periodic nodes if we are in a neighbouring tree
1825  bool is_periodic=false;
1826  if(in_neighbouring_tree)
1827  {
1828  //Is it periodic
1829  is_periodic =
1830  tree_pt()->root_pt()->is_neighbour_periodic(edges[edge_counter]);
1831  }
1832 
1833  // Loop over nodes along the edge
1834  for(unsigned i0=0;i0<n_p;i0++)
1835  {
1836  //Storage for pointer to the local node
1837  Node* local_node_pt=0;
1838 
1839  switch(edge_counter)
1840  {
1841  case 0:
1842  // Local fraction of node
1843  s_fraction[0] = local_one_d_fraction_of_node(i0,0);
1844  s_fraction[1] = 0.0;
1845  // Get pointer to local node
1846  local_node_pt = node_pt(i0);
1847  break;
1848 
1849  case 1:
1850  // Local fraction of node
1851  s_fraction[0] = local_one_d_fraction_of_node(i0,0);
1852  s_fraction[1] = 1.0;
1853  // Get pointer to local node
1854  local_node_pt = node_pt(i0 + n_p*(n_p-1));
1855  break;
1856 
1857  case 2:
1858  // Local fraction of node
1859  s_fraction[0] = 0.0;
1860  s_fraction[1] = local_one_d_fraction_of_node(i0,1);
1861  // Get pointer to local node
1862  local_node_pt = node_pt(n_p*i0);
1863  break;
1864 
1865  case 3:
1866  // Local fraction of node
1867  s_fraction[0] = 1.0;
1868  s_fraction[1] = local_one_d_fraction_of_node(i0,1);
1869  // Get pointer to local node
1870  local_node_pt = node_pt(n_p-1 + n_p*i0);
1871  break;
1872  }
1873 
1874  //Calculate the local coordinate and the local coordinate as viewed
1875  //from the neighbour
1876  Vector<double> s_in_neighb(2);
1877  for(unsigned i=0;i<2;i++)
1878  {
1879  //Local coordinate in this element
1880  s[i] = -1.0 + 2.0*s_fraction[i];
1881  //Local coordinate in the neighbour
1882  s_in_neighb[i] = s_lo_neigh[i] + s_fraction[translate_s[i]]*
1883  (s_hi_neigh[i] - s_lo_neigh[i]);
1884  }
1885 
1886  //Loop over timesteps
1887  for(unsigned t=0;t<n_time;t++)
1888  {
1889  // Get the nodal position from neighbour element
1890  Vector<double> x_in_neighb(2);
1891  neigh_pt->object_pt()->interpolated_x(t,s_in_neighb,x_in_neighb);
1892 
1893  // Check error only if the node is NOT periodic
1894  if(is_periodic==false)
1895  {
1896  for(int i=0;i<2;i++)
1897  {
1898  //Find the spatial error
1899  double err = std::fabs(local_node_pt->x(t,i) - x_in_neighb[i]);
1900 
1901  //If it's bigger than our tolerance, say so
1902  if (err>1e-9)
1903  {
1904  oomph_info << "errx " << err << " " << t << " "
1905  << local_node_pt->x(t,i)
1906  << " " << x_in_neighb[i]<< std::endl;
1907 
1908  oomph_info << "at " << local_node_pt->x(0) << " "
1909  << local_node_pt->x(1) << std::endl;
1910  }
1911 
1912  //If it's bigger than the previous max error, it is the
1913  //new max error!
1914  if (err>max_error_x[i]) {max_error_x[i]=err;}
1915  }
1916  }
1917 
1918  // Get the values from neighbour element. Note: # of values
1919  // gets set by routine (because in general we don't know
1920  // how many interpolated values a certain element has
1921  Vector<double> values_in_neighb;
1922  neigh_pt->object_pt()->
1923  get_interpolated_values(t,s_in_neighb,values_in_neighb);
1924 
1925  // Get the values in current element.
1926  Vector<double> values;
1927  get_interpolated_values(t,s,values);
1928 
1929  // Now figure out how many continuously interpolated values there are
1930  unsigned num_val=neigh_pt->object_pt()->ncont_interpolated_values();
1931 
1932  // Check error
1933  for(unsigned ival=0;ival<num_val;ival++)
1934  {
1935  double err=std::fabs(values[ival] - values_in_neighb[ival]);
1936 
1937  if (err>1.0e-10)
1938  {
1939  oomph_info << local_node_pt->x(0) << " "
1940  << local_node_pt->x(1) << " \n# "
1941  << "erru (S)" << err << " " << ival << " "
1942  << get_node_number(local_node_pt) << " "
1943  << values[ival]
1944  << " " << values_in_neighb[ival] << std::endl;
1945  }
1946 
1947  if (err>max_error_val) {max_error_val=err;}
1948 
1949  }
1950  }
1951 
1952  }
1953  }
1954  }
1955 
1956  max_error=max_error_x[0];
1957  if (max_error_x[1]>max_error) max_error=max_error_x[1];
1958  if (max_error_val>max_error) max_error=max_error_val;
1959 
1960  if (max_error>1e-9)
1961  {
1962  oomph_info << "\n#------------------------------------ \n#Max error " ;
1963  oomph_info << max_error_x[0]
1964  << " " << max_error_x[1]
1965  << " " << max_error_val << std::endl;
1966  oomph_info << "#------------------------------------ \n " << std::endl;
1967 
1968  }
1969 
1970 }
1971 
1972 
1973 
1974 //========================================================================
1975 /// Static matrix for coincidence between son nodal points and
1976 /// father boundaries
1977 ///
1978 //========================================================================
1979 std::map<unsigned,DenseMatrix<int> > RefineableQElement<2>::Father_bound;
1980 
1981 
1982 
1983 
1984 
1985 
1986 //////////////////////////////////////////////////////////////////////////
1987 //////////////////////////////////////////////////////////////////////////
1988 //////////////////////////////////////////////////////////////////////////
1989 
1990 
1991 
1992 
1993 //==================================================================
1994 /// Determine vector of solid (positional) boundary conditions along
1995 /// the element's boundary (or vertex) bound (S/W/N/E/SW/SE/NW/NE).
1996 ///
1997 /// This function assumes that the same boundary condition is applied
1998 /// along the entire length of an element's edge (of course, the
1999 /// vertices combine the boundary conditions of their two adjacent edges
2000 /// in the most restrictive combination. Hence, if we're at a vertex,
2001 /// we apply the most restrictive boundary condition of the
2002 /// two adjacent edges. If we're on an edge (in its proper interior),
2003 /// we apply the least restrictive boundary condition of all nodes
2004 /// along the edge.
2005 ///
2006 /// Usual convention:
2007 /// - solid_bound_cons[i]=0 if displacement in coordinate direction i
2008 /// on this boundary is free.
2009 /// - solid_bound_cons[i]=1 if it's pinned.
2010 //==================================================================
2012  Vector<int>& solid_bound_cons) const
2013 {
2014 
2015  using namespace QuadTreeNames;
2016 
2017  // Spatial dimension of all nodes
2018  unsigned n_dim = this->nodal_dimension();
2019 
2020  //Set up temporary vectors to hold edge boundary conditions
2021  Vector<int> bound_cons1(n_dim), bound_cons2(n_dim);
2022 
2023  //Which boundary are we on?
2024  switch(bound)
2025  {
2026  //If on edge, just get the bcs
2027  case N:
2028  case S:
2029  case W:
2030  case E:
2031  get_edge_solid_bcs(bound,solid_bound_cons);
2032  break;
2033 
2034  //Most restrictive boundary at SE corner
2035  case SE:
2036  get_edge_solid_bcs(S,bound_cons1);
2037  get_edge_solid_bcs(E,bound_cons2);
2038 
2039  for(unsigned k=0;k<n_dim;k++)
2040  {solid_bound_cons[k]=(bound_cons1[k]||bound_cons2[k]);}
2041  break;
2042 
2043  //Most restrictive boundary at SW corner
2044  case SW:
2045  get_edge_solid_bcs(S,bound_cons1);
2046  get_edge_solid_bcs(W,bound_cons2);
2047 
2048  for(unsigned k=0;k<n_dim;k++)
2049  {solid_bound_cons[k]=(bound_cons1[k]||bound_cons2[k]);}
2050  break;
2051 
2052  //Most restrictive boundary at NW corner
2053  case NW:
2054  get_edge_solid_bcs(N,bound_cons1);
2055  get_edge_solid_bcs(W,bound_cons2);
2056 
2057  for(unsigned k=0;k<n_dim;k++)
2058  {solid_bound_cons[k]=(bound_cons1[k]||bound_cons2[k]);}
2059  break;
2060 
2061  //Most restrictive boundary at NE corner
2062  case NE:
2063  get_edge_solid_bcs(N,bound_cons1);
2064  get_edge_solid_bcs(E,bound_cons2);
2065 
2066  for(unsigned k=0;k<n_dim;k++)
2067  {solid_bound_cons[k]=(bound_cons1[k]||bound_cons2[k]);}
2068  break;
2069 
2070  default:
2071  throw OomphLibError("Wrong boundary",
2072  OOMPH_CURRENT_FUNCTION,
2073  OOMPH_EXCEPTION_LOCATION);
2074  }
2075 }
2076 
2077 //==================================================================
2078 /// Determine Vector of solid (positional) boundary conditions along
2079 /// the element's edge (S/N/W/E) -- BC is the least restrictive combination
2080 /// of all the nodes on this edge
2081 ///
2082 /// Usual convention:
2083 /// - solid_bound_cons[i]=0 if displacement in coordinate direction i
2084 /// on this boundary is free
2085 /// - solid_bound_cons[i]=1 if it's pinned
2086 //==================================================================
2088  Vector<int>& solid_bound_cons) const
2089 {
2090  using namespace QuadTreeNames;
2091 
2092  //Number of nodes along 1D edge
2093  unsigned n_p = nnode_1d();
2094 
2095  //Left- and Right-hand nodes
2096  unsigned left_node, right_node;
2097 
2098  //Set the left (lower) and right (upper) nodes for the edge
2099  switch(edge)
2100  {
2101  case N:
2102  left_node = n_p*(n_p-1);
2103  right_node = n_p*n_p - 1;
2104  break;
2105 
2106  case S:
2107  left_node = 0;
2108  right_node = n_p-1;
2109  break;
2110 
2111  case W:
2112  left_node = 0;
2113  right_node = n_p*(n_p-1);
2114  break;
2115 
2116  case E:
2117  left_node = n_p-1;
2118  right_node = n_p*n_p - 1;
2119  break;
2120 
2121  default:
2122  std::ostringstream error_stream;
2123  error_stream
2124  << "Wrong edge " << edge << " passed to get_solid_edge_bcs(..)"
2125  << std::endl;
2126 
2127  throw OomphLibError(error_stream.str(),
2128  OOMPH_CURRENT_FUNCTION,
2129  OOMPH_EXCEPTION_LOCATION);
2130  }
2131 
2132 
2133  // Cast to SolidNodes
2134  SolidNode* left_node_pt =dynamic_cast<SolidNode*>(node_pt(left_node));
2135  SolidNode* right_node_pt=dynamic_cast<SolidNode*>(node_pt(right_node));
2136 #ifdef PARANOID
2137  if (left_node_pt==0)
2138  {
2139  throw OomphLibError(
2140  "Left node cannot be cast to SolidNode --> something is wrong",
2141  OOMPH_CURRENT_FUNCTION,
2142  OOMPH_EXCEPTION_LOCATION);
2143  }
2144  if (right_node_pt==0)
2145  {
2146  throw OomphLibError(
2147  "Right node cannot be cast to SolidNode --> something is wrong",
2148  OOMPH_CURRENT_FUNCTION,
2149  OOMPH_EXCEPTION_LOCATION);
2150  }
2151 #endif
2152 
2153 
2154  // Number of coordinate directions
2155  unsigned n_dim = this->nodal_dimension();
2156 
2157  //Loop over all directions, the least restrictive value is
2158  //the boundary condition at the left multiplied by that at the right
2159  //Assuming that free is always zero and pinned is one
2160  for(unsigned k=0;k<n_dim;k++)
2161  {
2162  solid_bound_cons[k] =
2163  left_node_pt ->position_is_pinned(k)*
2164  right_node_pt->position_is_pinned(k);
2165  }
2166 }
2167 
2168 /// Build required templates
2169 //template class RefineableQElement<2>;
2170 
2171 }
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2301
virtual Node * interpolating_node_pt(const unsigned &n, const int &value_id)
In mixed elements, different sets of nodes are used to interpolate different unknowns. This function returns the n-th node that interpolates the value_id-th unknown. Default implementation is that all variables use the positional nodes, i.e. isoparametric elements. Note that any overloaded versions of this function MUST provide a set of nodes for the position, which always has the value_id -1.
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
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:246
cstr elem_len * i
Definition: cfortran.h:607
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get all continously interpolated function values in this element as a Vector. Note: Vector sets is ow...
bool is_neighbour_periodic(const int &direction)
Return whether the neighbour in the particular direction is periodic.
Definition: tree.h:350
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...
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
OomphInfo oomph_info
void set_master_node_pt(const unsigned &i, Node *const &master_node_pt, const double &weight)
Set the pointer to the i-th master node and its weight.
Definition: nodes.cc:1417
e
Definition: cfortran.h:575
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
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
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1693
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
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
virtual unsigned ninterpolating_node(const int &value_id)
Return the number of nodes that are used to interpolate the value_id-th unknown. Default is to assume...
void start(const unsigned &i)
(Re-)start i-th timer
void get_boundaries(const int &edge, std::set< unsigned > &boundaries) const
Determine set of (mesh) boundaries that the element edge/vertex lives on.
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
bool boundary_coordinate_exists(const unsigned &i) const
Indicate whether the i-th boundary has an intrinsic coordinate.
Definition: mesh.h:565
QuadTree * gteq_edge_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level, bool &in_neighbouring_tree) const
Return pointer to greater or equal-sized edge neighbour in specified direction; also provide info reg...
Definition: quadtree.cc:423
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.
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:549
Class that contains data for hanging nodes.
Definition: nodes.h:684
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
void set_hanging_pt(HangInfo *const &hang_pt, const int &i)
Set the hanging data for the i-th value. (hang_pt=0 to make non-hanging)
Definition: nodes.cc:1997
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
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
Base class for elements that allow MacroElement-based node update.
int son_type() const
Return son type.
Definition: tree.h:209
virtual void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
Return the basis functions that are used to interpolate the value_id-th unknown. By default assume is...
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.
Refineable version of Solid quad elements.
void get_bcs(int bound, Vector< int > &bound_cons) const
Determine Vector of boundary conditions along edge (or on vertex) bound (S/W/N/E/SW/SE/NW/NE): For va...
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
Vector< std::string > colour
Tecplot colours.
void interpolated_zeta_on_edge(const unsigned &boundary, const int &edge, const Vector< double > &s, Vector< double > &zeta)
Return the value of the intrinsic boundary coordinate interpolated along the edge (S/W/N/E) ...
int & method_for_shape_derivs()
Access to method (enumerated flag) for determination of shape derivs.
virtual double s_min() const
Min value of local coordinate.
Definition: elements.h:2631
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
TreeRoot *& root_pt()
Return pointer to root of the tree.
Definition: tree.h:143
MacroElement * Macro_elem_pt
Pointer to the element's macro element (NULL by default)
Definition: elements.h:1644