refineable_brick_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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #include <algorithm>
27 
28 #include "mesh.h"
29 #include "algebraic_elements.h"
30 #include "macro_element_node_update_element.h"
31 #include "refineable_brick_element.h"
32 
33 
34 namespace oomph
35 {
36  //========================================================================
37  /// Print corner nodes, use colour (default "BLACK")
38  /// in right order so that tecplot can draw a cube without crossed lines
39  //========================================================================
40  void RefineableQElement<3>::output_corners(std::ostream& outfile,
41  const std::string& colour) const
42  {
43  Vector<double> s(3);
44  Vector<double> corner(3);
45 
46  outfile << "ZONE I=2, J=2, K=2 C=" << colour << std::endl;
47 
48  s[0] = -1.0;
49  s[1] = -1.0;
50  s[2] = -1.0;
51  get_x(s, corner);
52  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
53  << Number << std::endl;
54 
55  s[0] = -1.0;
56  s[1] = -1.0;
57  s[2] = 1.0;
58  get_x(s, corner);
59  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
60  << Number << std::endl;
61 
62  s[0] = -1.0;
63  s[1] = 1.0;
64  s[2] = -1.0;
65  get_x(s, corner);
66  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
67  << Number << std::endl;
68 
69  s[0] = -1.0;
70  s[1] = 1.0;
71  s[2] = 1.0;
72  get_x(s, corner);
73  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
74  << Number << std::endl;
75 
76  // next face
77 
78 
79  s[0] = 1.0;
80  s[1] = -1.0;
81  s[2] = -1.0;
82  get_x(s, corner);
83  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
84  << Number << std::endl;
85 
86  s[0] = 1.0;
87  s[1] = -1.0;
88  s[2] = 1.0;
89  get_x(s, corner);
90  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
91  << Number << std::endl;
92 
93  s[0] = 1.0;
94  s[1] = 1.0;
95  s[2] = -1.0;
96  get_x(s, corner);
97  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
98  << Number << std::endl;
99 
100  s[0] = 1.0;
101  s[1] = 1.0;
102  s[2] = 1.0;
103  get_x(s, corner);
104  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
105  << Number << std::endl;
106 
107 
108  // outfile << "TEXT CS = GRID, X = " << corner[0] <<
109  // ",Y = " << corner[1] << ",Z = " << corner[2] <<
110  // ", HU = GRID, H = 0.01, AN = MIDCENTER, T =\""
111  // << Number << "\"" << std::endl;
112  }
113 
114 
115  //==================================================================
116  /// Setup static matrix for coincidence between son nodal points and
117  /// father boundaries:
118  ///
119  /// Father_bound[nnode_1d](nnode_son,son_type)={RU/RF/RD/RB/.../ OMEGA}
120  ///
121  /// so that node nnode_son in element of type son_type lies
122  /// on boundary/vertex Father_bound[nnode_1d](nnode_son,son_type) in its
123  /// father element. If the node doesn't lie on a boundary
124  /// the value is OMEGA.
125  //==================================================================
126  void RefineableQElement<3>::setup_father_bounds()
127  {
128  using namespace OcTreeNames;
129 
130  // Find the number of nodes along a 1D edge
131  unsigned n_p = nnode_1d();
132 
133  // Allocate space for the boundary information
134  Father_bound[n_p].resize(n_p * n_p * n_p, 8);
135 
136  // Initialise: By default points are not on the boundary
137  for (unsigned n = 0; n < n_p * n_p * n_p; n++)
138  {
139  for (unsigned ison = 0; ison < 8; ison++)
140  {
141  Father_bound[n_p](n, ison) = Tree::OMEGA;
142  }
143  }
144 
145  for (int i_son = LDB; i_son <= RUF; i_son++)
146  {
147  // vector representing the son
148  Vector<int> vect_son(3);
149  // vector representing (at the end) the boundaries
150  Vector<int> vect_bound(3);
151  vect_son = OcTree::Direction_to_vector[i_son];
152  for (unsigned i0 = 0; i0 < n_p; i0++)
153  {
154  for (unsigned i1 = 0; i1 < n_p; i1++)
155  {
156  for (unsigned i2 = 0; i2 < n_p; i2++)
157  {
158  // Initialisation to make it work
159  for (unsigned i = 0; i < 3; i++)
160  {
161  vect_bound[i] = -vect_son[i];
162  }
163  // Seting up the boundaries coordinates as if the coordinates
164  // of vect_bound had been initialised to 0 if it were so,
165  // vect_bound would be the vector of the boundaries in the son
166  // itself.
167 
168  if (i0 == 0)
169  {
170  vect_bound[0] = -1;
171  }
172  if (i0 == n_p - 1)
173  {
174  vect_bound[0] = 1;
175  }
176  if (i1 == 0)
177  {
178  vect_bound[1] = -1;
179  }
180  if (i1 == n_p - 1)
181  {
182  vect_bound[1] = 1;
183  }
184  if (i2 == 0)
185  {
186  vect_bound[2] = -1;
187  }
188  if (i2 == n_p - 1)
189  {
190  vect_bound[2] = 1;
191  }
192 
193  // The effect of this is to filter the boundaries to keep only the
194  // ones which are effectively father boundaries.
195  // -- if the node is not on a "i0 boundary", we still
196  // have vect_bound[0]=-vect_son[0]
197  // and the result is vect_bound[0]=0
198  // (he is not on this boundary for his father)
199  // -- if he is on a son's boundary which is not one of
200  // the father -> same thing
201  // -- if he is on a boundary which is one of his father's,
202  // vect_bound[i]=vect_son[i]
203  // and the new vect_bound[i] is the same as the old one
204  for (int i = 0; i < 3; i++)
205  {
206  vect_bound[i] = (vect_bound[i] + vect_son[i]) / 2;
207  }
208 
209  // Return the result as {U,R,D,...RDB,LUF,OMEGA}
210  Father_bound[n_p](i0 + n_p * i1 + n_p * n_p * i2, i_son) =
211  OcTree::Vector_to_direction[vect_bound];
212 
213  } // Loop over i2
214  } // Loop over i1
215  } // Loop over i0
216  } // Loop over i_son
217 
218  } // setup_father_bounds()
219 
220 
221  //==================================================================
222  /// Determine Vector of boundary conditions along the element's boundary
223  /// bound.
224  ///
225  /// This function assumes that the same boundary condition is applied
226  /// along the entire area of an element's face (of course, the
227  /// vertices combine the boundary conditions of their two adjacent edges
228  /// in the most restrictive combination. Hence, if we're at a vertex,
229  /// we apply the most restrictive boundary condition of the
230  /// two adjacent edges. If we're on an edge (in its proper interior),
231  /// we apply the least restrictive boundary condition of all nodes
232  /// along the edge.
233  ///
234  /// Usual convention:
235  /// - bound_cons[ival]=0 if value ival on this boundary is free
236  /// - bound_cons[ival]=1 if value ival on this boundary is pinned
237  //==================================================================
238  void RefineableQElement<3>::get_bcs(int bound, Vector<int>& bound_cons) const
239  {
240  using namespace OcTreeNames;
241 
242  // Max. number of nodal data values in element
243  unsigned nvalue = ncont_interpolated_values();
244  // Set up temporary vectors to hold edge boundary conditions
245  Vector<int> bound_cons1(nvalue), bound_cons2(nvalue);
246  Vector<int> bound_cons3(nvalue);
247 
248  Vector<int> vect1(3), vect2(3), vect3(3);
249  Vector<int> vect_elem;
250  Vector<int> notzero;
251  int n = 0;
252 
253  vect_elem = OcTree::Direction_to_vector[bound];
254 
255  // Just to see if bound is a face, an edge, or a vertex, n stores
256  // the number of non-zero values in the vector reprensenting the bound
257  // and the vector notzero stores the position of these values
258  for (int i = 0; i < 3; i++)
259  {
260  if (vect_elem[i] != 0)
261  {
262  n++;
263  notzero.push_back(i);
264  }
265  }
266 
267  switch (n)
268  {
269  // If there is only one non-zero value, bound is a face
270  case 1:
271  get_face_bcs(bound, bound_cons);
272  break;
273 
274  // If there are two non-zero values, bound is an edge
275  case 2:
276 
277  for (unsigned i = 0; i < 3; i++)
278  {
279  vect1[i] = 0;
280  vect2[i] = 0;
281  }
282  // vect1 and vect2 are the vector of the two faces adjacent to bound
283  vect1[notzero[0]] = vect_elem[notzero[0]];
284  vect2[notzero[1]] = vect_elem[notzero[1]];
285 
286  get_face_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
287  get_face_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
288  // get the most restrictive bc
289  for (unsigned k = 0; k < nvalue; k++)
290  {
291  bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
292  }
293  break;
294 
295  // If there are three non-zero value, bound is a vertex
296  case 3:
297 
298  for (unsigned i = 0; i < 3; i++)
299  {
300  vect1[i] = 0;
301  vect2[i] = 0;
302  vect3[i] = 0;
303  }
304  // vectors to the three adjacent faces of the vertex
305  vect1[0] = vect_elem[0];
306  vect2[1] = vect_elem[1];
307  vect3[2] = vect_elem[2];
308 
309  get_face_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
310  get_face_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
311  get_face_bcs(OcTree::Vector_to_direction[vect3], bound_cons3);
312 
313 
314  // set the bcs to the most restrictive ones
315  for (unsigned k = 0; k < nvalue; k++)
316  {
317  bound_cons[k] = (bound_cons1[k] || bound_cons2[k] || bound_cons3[k]);
318  }
319  break;
320 
321  default:
322  throw OomphLibError("Make sure you are not giving OMEGA as bound",
323  OOMPH_CURRENT_FUNCTION,
324  OOMPH_EXCEPTION_LOCATION);
325  }
326  }
327 
328  //==================================================================
329  /// Determine Vector of boundary conditions along the element's
330  /// face (R/L/U/D/B/F) -- BC is the least restrictive combination
331  /// of all the nodes on this face
332  ///
333  /// Usual convention:
334  /// - bound_cons[ival]=0 if value ival on this boundary is free
335  /// - bound_cons[ival]=1 if value ival on this boundary is pinned
336  //==================================================================
337  void RefineableQElement<3>::get_face_bcs(const int& face,
338  Vector<int>& bound_cons) const
339  {
340  using namespace OcTreeNames;
341 
342  // Number of nodes along 1D edge
343  unsigned n_p = nnode_1d();
344  // the four corner nodes on the boundary
345  unsigned node1, node2, node3, node4;
346 
347  // Set the four corner nodes for the face
348  switch (face)
349  {
350  case U:
351  node1 = n_p * n_p * n_p - 1;
352  node2 = n_p * n_p - 1;
353  node3 = n_p * (n_p - 1);
354  node4 = n_p * (n_p * n_p - 1);
355  break;
356 
357  case D:
358  node1 = 0;
359  node2 = n_p - 1;
360  node3 = (n_p * n_p + 1) * (n_p - 1);
361  node4 = n_p * n_p * (n_p - 1);
362  break;
363 
364  case R:
365  node1 = n_p - 1;
366  node2 = (n_p * n_p + 1) * (n_p - 1);
367  node3 = n_p * n_p * n_p - 1;
368  node4 = n_p * n_p - 1;
369  break;
370 
371  case L:
372  node1 = 0;
373  node2 = n_p * (n_p - 1);
374  node3 = n_p * (n_p * n_p - 1);
375  node4 = n_p * n_p * (n_p - 1);
376  break;
377 
378  case B:
379  node1 = 0;
380  node2 = n_p - 1;
381  node3 = n_p * n_p - 1;
382  node4 = n_p * (n_p - 1);
383  break;
384 
385  case F:
386  node1 = n_p * n_p * n_p - 1;
387  node2 = n_p * (n_p * n_p - 1);
388  node3 = n_p * n_p * (n_p - 1);
389  node4 = (n_p - 1) * (n_p * n_p + 1);
390  break;
391 
392  default:
393  std::ostringstream error_stream;
394  error_stream << "Wrong edge " << face << " passed\n";
395 
396  throw OomphLibError(
397  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
398  }
399 
400  // Max. number of nodal data values in element
401  unsigned maxnvalue = ncont_interpolated_values();
402 
403  // Loop over all values, the least restrictive value is
404  // the multiplication of the boundary conditions at the 4 nodes
405  // Assuming that free is always zero and pinned is one
406  for (unsigned k = 0; k < maxnvalue; k++)
407  {
408  bound_cons[k] =
409  node_pt(node1)->is_pinned(k) * node_pt(node2)->is_pinned(k) *
410  node_pt(node3)->is_pinned(k) * node_pt(node4)->is_pinned(k);
411  }
412  }
413 
414 
415  //==================================================================
416  /// Given an element edge/vertex, return a Vector which contains
417  /// all the (mesh-)boundary numbers that this element edge/vertex
418  /// lives on.
419  ///
420  /// For proper edges, the boundary is the one (if any) that is shared by
421  /// both vertex nodes). For vertex nodes, we just return their
422  /// boundaries.
423  //==================================================================
424  void RefineableQElement<3>::get_boundaries(const int& element,
425  std::set<unsigned>& boundary) const
426  {
427  using namespace OcTreeNames;
428 
429  // Number of 1d nodes along an edge
430  unsigned n_p = nnode_1d();
431  // Left and right-hand nodes
432  int node[4];
433  int a = 0, b = 0;
434  int n = 0;
435  Vector<int> vect_face(3);
436  vect_face = OcTree::Direction_to_vector[element];
437  // Set the left (lower) and right (upper) nodes for the edge
438 
439  // this is to know what is the type of element (face, edge, vertex)
440  // we just need to count the number of values equal to 0 in the
441  // vector representing this element
442 
443  // Local storage for the node-numbers in given directions
444  // Initialise to zero (assume at LH end of each domain)
445  int one_d_node_number[3] = {0, 0, 0};
446  // n stores the number of values equal to 0, a is the position of the
447  // last 0-value ;b is the position of the last non0-value
448  for (int i = 0; i < 3; i++)
449  {
450  if (vect_face[i] == 0)
451  {
452  a = i;
453  n++;
454  }
455  else
456  {
457  b = i;
458  // If we are at the extreme (RH) end of the face,
459  // set the node number accordingly
460  if (vect_face[i] == 1)
461  {
462  one_d_node_number[i] = n_p - 1;
463  }
464  }
465  }
466 
467  switch (n)
468  {
469  // if n=0 element is a vertex, and need to look at only one node
470  case 0:
471  node[0] = one_d_node_number[0] + n_p * one_d_node_number[1] +
472  n_p * n_p * one_d_node_number[2];
473  node[1] = node[0];
474  node[2] = node[0];
475  node[3] = node[0];
476  break;
477 
478  // if n=1 element is an edge, and we need to look at two nodes
479  case 1:
480  if (a == 0)
481  {
482  node[0] = (n_p - 1) + n_p * one_d_node_number[1] +
483  n_p * n_p * one_d_node_number[2];
484  node[1] =
485  n_p * one_d_node_number[1] + n_p * n_p * one_d_node_number[2];
486  }
487  else if (a == 1)
488  {
489  node[0] = n_p * (n_p - 1) + one_d_node_number[0] +
490  n_p * n_p * one_d_node_number[2];
491  node[1] = one_d_node_number[0] + n_p * n_p * one_d_node_number[2];
492  }
493  else if (a == 2)
494  {
495  node[0] = one_d_node_number[0] + n_p * one_d_node_number[1] +
496  n_p * n_p * (n_p - 1);
497  node[1] = one_d_node_number[0] + n_p * one_d_node_number[1];
498  }
499  node[2] = node[1];
500  node[3] = node[1];
501  break;
502 
503  // if n=2 element is a face, and we need to look at its 4 nodes
504  case 2:
505  if (b == 0)
506  {
507  node[0] =
508  one_d_node_number[0] + n_p * n_p * (n_p - 1) + n_p * (n_p - 1);
509  node[1] = one_d_node_number[0] + n_p * (n_p - 1);
510  node[2] = one_d_node_number[0] + n_p * n_p * (n_p - 1);
511  node[3] = one_d_node_number[0];
512  }
513  else if (b == 1)
514  {
515  node[0] =
516  n_p * one_d_node_number[1] + n_p * n_p * (n_p - 1) + (n_p - 1);
517  node[1] = n_p * one_d_node_number[1] + (n_p - 1);
518  node[2] = n_p * one_d_node_number[1] + n_p * n_p * (n_p - 1);
519  node[3] = n_p * one_d_node_number[1];
520  }
521  else if (b == 2)
522  {
523  node[0] =
524  n_p * n_p * one_d_node_number[2] + n_p * (n_p - 1) + (n_p - 1);
525  node[1] = n_p * n_p * one_d_node_number[2] + (n_p - 1);
526  node[2] = n_p * n_p * one_d_node_number[2] + n_p * (n_p - 1);
527  node[3] = n_p * n_p * one_d_node_number[2];
528  }
529  break;
530  default:
531  throw OomphLibError("Make sure you are not giving OMEGA as boundary",
532  OOMPH_CURRENT_FUNCTION,
533  OOMPH_EXCEPTION_LOCATION);
534  }
535 
536 
537  // Empty boundary set: Edge does not live on any boundary
538  boundary.clear();
539 
540  // Storage for the boundaries at the four nodes
541  Vector<std::set<unsigned>*> node_boundaries_pt(4, 0);
542 
543  // Loop over the four nodes and get the boundary information
544  for (unsigned i = 0; i < 4; i++)
545  {
546  node_pt(node[i])->get_boundaries_pt(node_boundaries_pt[i]);
547  }
548 
549 
550  // Now work out the intersections
551  Vector<std::set<unsigned>> boundary_aux(2);
552 
553  for (unsigned i = 0; i < 2; i++)
554  {
555  // If the two nodes both lie on boundaries
556  if ((node_boundaries_pt[2 * i] != 0) &&
557  (node_boundaries_pt[2 * i + 1] != 0))
558  {
559  // Find the intersection (the common boundaries) of these nodes
560  std::set_intersection(
561  node_boundaries_pt[2 * i]->begin(),
562  node_boundaries_pt[2 * i]->end(),
563  node_boundaries_pt[2 * i + 1]->begin(),
564  node_boundaries_pt[2 * i + 1]->end(),
565  inserter(boundary_aux[i], boundary_aux[i].begin()));
566  }
567  }
568 
569  // Now calculate the total intersection
570  set_intersection(boundary_aux[0].begin(),
571  boundary_aux[0].end(),
572  boundary_aux[1].begin(),
573  boundary_aux[1].end(),
574  inserter(boundary, boundary.begin()));
575  }
576 
577 
578  //===================================================================
579  /// Return the value of the intrinsic boundary coordinate interpolated
580  /// along the face
581  //===================================================================
582  void RefineableQElement<3>::interpolated_zeta_on_face(
583  const unsigned& boundary,
584  const int& face,
585  const Vector<double>& s,
586  Vector<double>& zeta)
587  {
588  using namespace OcTreeNames;
589 
590  // Number of nodes along an edge
591  unsigned nnodes_1d = nnode_1d();
592 
593  // Number of nodes on a face
594  unsigned nnodes_2d = nnodes_1d * nnodes_1d;
595 
596  // Total number of nodes
597  unsigned nnodes_3d = nnode();
598 
599  // Storage for the shape functions
600  Shape psi(nnodes_3d);
601 
602  // Get the shape functions at the passed position
603  this->shape(s, psi);
604 
605  // Unsigned data that give starts and increments for the loop
606  // over the nodes on the faces.
607  unsigned start = 0, increment1 = 1, increment2 = 1;
608 
609  // Flag to record if actually on a face or an edge
610  bool on_edge = true;
611 
612  // Which face?
613  switch (face)
614  {
615  case L:
616 #ifdef PARANOID
617  if (s[0] != -1.0)
618  {
619  std::ostringstream error_stream;
620  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
621  << " is not on Left face\n";
622 
623  throw OomphLibError(error_stream.str(),
624  OOMPH_CURRENT_FUNCTION,
625  OOMPH_EXCEPTION_LOCATION);
626  }
627 #endif
628  // Start is zero (bottom-left-back corner)
629  increment1 = nnodes_1d;
630  increment2 = 0;
631  on_edge = false;
632  break;
633 
634  case R:
635 #ifdef PARANOID
636  if (s[0] != 1.0)
637  {
638  std::ostringstream error_stream;
639  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
640  << " is not on Right face\n";
641 
642  throw OomphLibError(error_stream.str(),
643  OOMPH_CURRENT_FUNCTION,
644  OOMPH_EXCEPTION_LOCATION);
645  }
646 #endif
647  // Start is bottom-right-back corner
648  start = nnodes_1d - 1;
649  increment1 = nnodes_1d;
650  increment2 = 0;
651  on_edge = false;
652  break;
653 
654  case D:
655 #ifdef PARANOID
656  if (s[1] != -1.0)
657  {
658  std::ostringstream error_stream;
659  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
660  << " is not on Bottom face\n";
661 
662  throw OomphLibError(error_stream.str(),
663  OOMPH_CURRENT_FUNCTION,
664  OOMPH_EXCEPTION_LOCATION);
665  }
666 #endif
667  // Start is zero and increments2 is nnode_2d-nnode_1d
668  increment2 = nnodes_2d - nnodes_1d;
669  on_edge = false;
670  break;
671 
672  case U:
673 #ifdef PARANOID
674  if (s[1] != 1.0)
675  {
676  std::ostringstream error_stream;
677  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
678  << " is not on Upper face\n";
679 
680  throw OomphLibError(error_stream.str(),
681  OOMPH_CURRENT_FUNCTION,
682  OOMPH_EXCEPTION_LOCATION);
683  }
684 #endif
685  // Start is top-left-back corner and increments2 is nnode_2d-nnode_1d
686  start = nnodes_2d - nnodes_1d;
687  increment2 = start;
688  on_edge = false;
689  break;
690 
691  case B:
692 #ifdef PARANOID
693  if (s[2] != -1.0)
694  {
695  std::ostringstream error_stream;
696  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
697  << " is not on Back face\n";
698 
699  throw OomphLibError(error_stream.str(),
700  OOMPH_CURRENT_FUNCTION,
701  OOMPH_EXCEPTION_LOCATION);
702  }
703 #endif
704  // Start is zero and increments are 1 and 0
705  increment2 = 0;
706  on_edge = false;
707  break;
708 
709  case F:
710 #ifdef PARANOID
711  if (s[2] != 1.0)
712  {
713  std::ostringstream error_stream;
714  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
715  << " is not on Front face\n";
716 
717  throw OomphLibError(error_stream.str(),
718  OOMPH_CURRENT_FUNCTION,
719  OOMPH_EXCEPTION_LOCATION);
720  }
721 #endif
722  // Start is bottom-left-front corner
723  start = nnodes_3d - nnodes_2d;
724  increment2 = 0;
725  on_edge = false;
726  break;
727 
728  case LF:
729 #ifdef PARANOID
730  if ((s[0] != -1.0) || (s[2] != 1.0))
731  {
732  std::ostringstream error_stream;
733  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
734  << " is not on Front-Left edge\n";
735 
736  throw OomphLibError(error_stream.str(),
737  OOMPH_CURRENT_FUNCTION,
738  OOMPH_EXCEPTION_LOCATION);
739  }
740 #endif
741  start = nnodes_3d - nnodes_2d;
742  increment1 = nnodes_1d;
743  break;
744 
745  case LD:
746 #ifdef PARANOID
747  if ((s[0] != -1.0) || (s[1] != -1.0))
748  {
749  std::ostringstream error_stream;
750  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
751  << " is not on Bottom-Left edge\n";
752 
753  throw OomphLibError(error_stream.str(),
754  OOMPH_CURRENT_FUNCTION,
755  OOMPH_EXCEPTION_LOCATION);
756  }
757 #endif
758  increment1 = nnodes_2d;
759  break;
760 
761  case LU:
762 #ifdef PARANOID
763  if ((s[0] != -1.0) || (s[1] != 1.0))
764  {
765  std::ostringstream error_stream;
766  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
767  << " is not on Upper-Left edge\n";
768 
769  throw OomphLibError(error_stream.str(),
770  OOMPH_CURRENT_FUNCTION,
771  OOMPH_EXCEPTION_LOCATION);
772  }
773 #endif
774  start = nnodes_2d - nnodes_1d;
775  increment1 = nnodes_2d;
776  break;
777 
778  case LB:
779 #ifdef PARANOID
780  if ((s[0] != -1.0) || (s[2] != -1.0))
781  {
782  std::ostringstream error_stream;
783  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
784  << " is not on Back-Left edge\n";
785 
786  throw OomphLibError(error_stream.str(),
787  OOMPH_CURRENT_FUNCTION,
788  OOMPH_EXCEPTION_LOCATION);
789  }
790 #endif
791  increment1 = nnodes_1d;
792  break;
793 
794  case RF:
795 #ifdef PARANOID
796  if ((s[0] != 1.0) || (s[2] != 1.0))
797  {
798  std::ostringstream error_stream;
799  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
800  << " is not on Front-Right edge\n";
801 
802  throw OomphLibError(error_stream.str(),
803  OOMPH_CURRENT_FUNCTION,
804  OOMPH_EXCEPTION_LOCATION);
805  }
806 #endif
807  start = nnodes_3d - 1;
808  increment1 = -nnodes_1d;
809  break;
810 
811  case RD:
812 #ifdef PARANOID
813  if ((s[0] != 1.0) || (s[1] != -1.0))
814  {
815  std::ostringstream error_stream;
816  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
817  << " is not on Bottom-Right edge\n";
818 
819  throw OomphLibError(error_stream.str(),
820  OOMPH_CURRENT_FUNCTION,
821  OOMPH_EXCEPTION_LOCATION);
822  }
823 #endif
824  start = nnodes_1d - 1;
825  increment1 = nnodes_2d;
826  break;
827 
828  case RU:
829 #ifdef PARANOID
830  if ((s[0] != 1.0) || (s[1] != 1.0))
831  {
832  std::ostringstream error_stream;
833  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
834  << " is not on Upper-Right edge\n";
835 
836  throw OomphLibError(error_stream.str(),
837  OOMPH_CURRENT_FUNCTION,
838  OOMPH_EXCEPTION_LOCATION);
839  }
840 #endif
841  start = nnodes_2d - 1;
842  increment1 = nnodes_2d;
843  break;
844 
845  case RB:
846 #ifdef PARANOID
847  if ((s[0] != 1.0) || (s[2] != -1.0))
848  {
849  std::ostringstream error_stream;
850  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
851  << " is not on Back-Right edge\n";
852 
853  throw OomphLibError(error_stream.str(),
854  OOMPH_CURRENT_FUNCTION,
855  OOMPH_EXCEPTION_LOCATION);
856  }
857 #endif
858  start = nnodes_1d - 1;
859  increment1 = nnodes_1d;
860  break;
861 
862  case DB:
863 #ifdef PARANOID
864  if ((s[1] != -1.0) || (s[2] != -1.0))
865  {
866  std::ostringstream error_stream;
867  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
868  << " is not on Back-Bottom edge\n";
869 
870  throw OomphLibError(error_stream.str(),
871  OOMPH_CURRENT_FUNCTION,
872  OOMPH_EXCEPTION_LOCATION);
873  }
874 #endif
875  break;
876 
877  case DF:
878 #ifdef PARANOID
879  if ((s[1] != -1.0) || (s[2] != 1.0))
880  {
881  std::ostringstream error_stream;
882  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
883  << " is not on Front-Bottom edge\n";
884 
885  throw OomphLibError(error_stream.str(),
886  OOMPH_CURRENT_FUNCTION,
887  OOMPH_EXCEPTION_LOCATION);
888  }
889 #endif
890  start = nnodes_3d - nnodes_2d;
891  break;
892 
893  case UB:
894 #ifdef PARANOID
895  if ((s[1] != 1.0) || (s[2] != -1.0))
896  {
897  std::ostringstream error_stream;
898  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
899  << " is not on Back-Upper edge\n";
900 
901  throw OomphLibError(error_stream.str(),
902  OOMPH_CURRENT_FUNCTION,
903  OOMPH_EXCEPTION_LOCATION);
904  }
905 #endif
906  start = nnodes_2d - nnodes_1d;
907  break;
908 
909  case UF:
910 #ifdef PARANOID
911  if ((s[1] != 1.0) || (s[2] != 1.0))
912  {
913  std::ostringstream error_stream;
914  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
915  << " is not on Upper-Front edge\n";
916 
917  throw OomphLibError(error_stream.str(),
918  OOMPH_CURRENT_FUNCTION,
919  OOMPH_EXCEPTION_LOCATION);
920  }
921 #endif
922  start = nnodes_3d - nnodes_1d;
923  break;
924 
925  default:
926 
927  std::ostringstream error_stream;
928  error_stream << "Wrong face " << OcTree::Direct_string[face]
929  << " passed" << std::endl;
930  error_stream << "Trouble at : s= [" << s[0] << " " << s[1] << " "
931  << s[2] << "]\n";
932  Vector<double> x(3);
933  this->interpolated_x(s, x);
934  error_stream << "corresponding to : x= [" << x[0] << " " << x[1] << " "
935  << x[2] << "]\n";
936  throw OomphLibError(
937  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
938  }
939 
940  // Initialise the intrinsic coordinate
941  zeta[0] = 0.0;
942  zeta[1] = 0.0;
943 
944  // Set the start node number
945  unsigned node = start;
946 
947  if (on_edge)
948  {
949  for (unsigned l1 = 0; l1 < nnodes_1d; l1++)
950  {
951  // Get the intrinsic coordinate
952  Vector<double> node_zeta(2);
953  node_pt(node)->get_coordinates_on_boundary(boundary, node_zeta);
954 
955  // Now multiply by the shape function
956  zeta[0] += node_zeta[0] * psi(node);
957  zeta[1] += node_zeta[1] * psi(node);
958 
959  // Update node
960  node += increment1;
961  }
962  }
963  else
964  {
965  for (unsigned l2 = 0; l2 < nnodes_1d; l2++)
966  {
967  for (unsigned l1 = 0; l1 < nnodes_1d; l1++)
968  {
969  // Get the intrinsic coordinate
970  Vector<double> node_zeta(2);
971  node_pt(node)->get_coordinates_on_boundary(boundary, node_zeta);
972 
973  // Now multiply by the shape function
974  zeta[0] += node_zeta[0] * psi(node);
975  zeta[1] += node_zeta[1] * psi(node);
976 
977  // Update node
978  node += increment1;
979  }
980  // Update node
981  node += increment2;
982  }
983  }
984  }
985 
986  //===================================================================
987  /// If a neighbouring element has already created a node at a
988  /// position corresponding to the local fractional position within the
989  /// present element, s_fraction, return
990  /// a pointer to that node. If not, return NULL (0).
991  //===================================================================
992  Node* RefineableQElement<3>::node_created_by_neighbour(
993  const Vector<double>& s_fraction, bool& is_periodic)
994  {
995  using namespace OcTreeNames;
996 
997  // Calculate the faces/edges on which the node lies
998  Vector<int> faces;
999  Vector<int> edges;
1000 
1001  if (s_fraction[0] == 0.0)
1002  {
1003  faces.push_back(L);
1004  if (s_fraction[1] == 0.0)
1005  {
1006  edges.push_back(LD);
1007  }
1008  if (s_fraction[2] == 0.0)
1009  {
1010  edges.push_back(LB);
1011  }
1012  if (s_fraction[1] == 1.0)
1013  {
1014  edges.push_back(LU);
1015  }
1016  if (s_fraction[2] == 1.0)
1017  {
1018  edges.push_back(LF);
1019  }
1020  }
1021 
1022  if (s_fraction[0] == 1.0)
1023  {
1024  faces.push_back(R);
1025  if (s_fraction[1] == 0.0)
1026  {
1027  edges.push_back(RD);
1028  }
1029  if (s_fraction[2] == 0.0)
1030  {
1031  edges.push_back(RB);
1032  }
1033  if (s_fraction[1] == 1.0)
1034  {
1035  edges.push_back(RU);
1036  }
1037  if (s_fraction[2] == 1.0)
1038  {
1039  edges.push_back(RF);
1040  }
1041  }
1042 
1043  if (s_fraction[1] == 0.0)
1044  {
1045  faces.push_back(D);
1046  if (s_fraction[2] == 0.0)
1047  {
1048  edges.push_back(DB);
1049  }
1050  if (s_fraction[2] == 1.0)
1051  {
1052  edges.push_back(DF);
1053  }
1054  }
1055 
1056  if (s_fraction[1] == 1.0)
1057  {
1058  faces.push_back(U);
1059  if (s_fraction[2] == 0.0)
1060  {
1061  edges.push_back(UB);
1062  }
1063  if (s_fraction[2] == 1.0)
1064  {
1065  edges.push_back(UF);
1066  }
1067  }
1068 
1069  if (s_fraction[2] == 0.0)
1070  {
1071  faces.push_back(B);
1072  }
1073 
1074  if (s_fraction[2] == 1.0)
1075  {
1076  faces.push_back(F);
1077  }
1078 
1079  // Find the number of faces
1080  unsigned n_face = faces.size();
1081 
1082  // Find the number of edges
1083  unsigned n_edge = edges.size();
1084 
1085  Vector<unsigned> translate_s(3);
1086  Vector<double> s_lo_neigh(3);
1087  Vector<double> s_hi_neigh(3);
1088  Vector<double> s(3);
1089 
1090  int neigh_face, diff_level;
1091 
1092  // Loop over the faces on which the node lies
1093  //-------------------------------------------
1094  for (unsigned j = 0; j < n_face; j++)
1095  {
1096  // Boolean to indicate whether or not the neighbour has a
1097  // different Tree root
1098  bool in_neighbouring_tree;
1099 
1100  // Find pointer to neighbouring element along face
1101  OcTree* neigh_pt;
1102  neigh_pt = octree_pt()->gteq_face_neighbour(faces[j],
1103  translate_s,
1104  s_lo_neigh,
1105  s_hi_neigh,
1106  neigh_face,
1107  diff_level,
1108  in_neighbouring_tree);
1109 
1110  // Neighbour exists
1111  if (neigh_pt != 0)
1112  {
1113  // Have its nodes been created yet?
1114  if (neigh_pt->object_pt()->nodes_built())
1115  {
1116  // We now need to translate the nodal location, defined in terms
1117  // of the fractional coordinates of the present element into
1118  // those of its neighbour. For this we use the information returned
1119  // to use from the octree function.
1120 
1121  // Calculate the local coordinate in the neighbour
1122  // Note that we need to use the translation scheme in case
1123  // the local coordinates are swapped in the neighbour.
1124  for (unsigned i = 0; i < 3; i++)
1125  {
1126  s[i] = s_lo_neigh[i] +
1127  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
1128  }
1129 
1130  // Find the node in the neighbour
1131  Node* neighbour_node_pt =
1132  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
1133 
1134  // If there is a node, return it
1135  if (neighbour_node_pt != 0)
1136  {
1137  // Now work out whether it's a periodic boundary. This is
1138  // only possible if we have moved into a neighbouring tree
1139  if (in_neighbouring_tree)
1140  {
1141  // Return whether the neighbour is periodic
1142  is_periodic =
1143  octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
1144  }
1145 
1146  // Return the neighbour node pointer
1147  return neighbour_node_pt;
1148  }
1149  } // if (neigh_pt->object_pt()->nodes_built())
1150  } // if (neigh_pt!=0)
1151  } // for (unsigned j=0;j<n_face;j++)
1152 
1153  // Loop over the edges on which the node lies
1154  //------------------------------------------
1155  for (unsigned j = 0; j < n_edge; j++)
1156  {
1157  // Even if we restrict ourselves to true edge neighbours (i.e.
1158  // elements that are not also face neighbours) there may be multiple
1159  // edge neighbours across the edges between multiple root octrees.
1160  // When making the first call to OcTree::gteq_true_edge_neighbour(...)
1161  // we simply return the first one of these multiple edge neighbours
1162  // (if there are any at all, of course) and also return the total number
1163  // of true edge neighbours. If the node in question already exists
1164  // on the first edge neighbour we're done. If it doesn't it may exist
1165  // on other edge neighbours so we repeat the process over all
1166  // other edge neighbours (bailing out if a node is found, of course).
1167 
1168  // Initially return the zero-th true edge neighbour
1169  unsigned i_root_edge_neighbour = 0;
1170 
1171  // Initialise the total number of true edge neighbours
1172  unsigned nroot_edge_neighbour = 0;
1173 
1174  // Keep searching until we've found the node or until we've checked
1175  // all available edge neighbours
1176  bool keep_searching = true;
1177  while (keep_searching)
1178  {
1179  // Find pointer to neighbouring element along edge
1180  OcTree* neigh_pt;
1181  neigh_pt = octree_pt()->gteq_true_edge_neighbour(edges[j],
1182  i_root_edge_neighbour,
1183  nroot_edge_neighbour,
1184  translate_s,
1185  s_lo_neigh,
1186  s_hi_neigh,
1187  neigh_face,
1188  diff_level);
1189 
1190  // Neighbour exists
1191  if (neigh_pt != 0)
1192  {
1193  // Have its nodes been created yet?
1194  if (neigh_pt->object_pt()->nodes_built())
1195  {
1196  // We now need to translate the nodal location, defined in terms
1197  // of the fractional coordinates of the present element into
1198  // those of its neighbour. For this we use the information returned
1199  // to use from the octree function.
1200 
1201  // Calculate the local coordinate in the neighbour
1202  // Note that we need to use the translation scheme in case
1203  // the local coordinates are swapped in the neighbour.
1204  for (unsigned i = 0; i < 3; i++)
1205  {
1206  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
1207  (s_hi_neigh[i] - s_lo_neigh[i]);
1208  }
1209 
1210  // Find the node in the neighbour
1211  Node* neighbour_node_pt =
1212  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
1213 
1214  // If there is a node, return it
1215  if (neighbour_node_pt != 0)
1216  {
1217  // Get the faces on which the edge lies
1218  Vector<int> faces_attached_to_edge =
1219  OcTree::faces_of_common_edge(edges[j]);
1220 
1221  // Get the number of entries in the vector
1222  unsigned n_faces_attached_to_edge = faces_attached_to_edge.size();
1223 
1224  // Loop over the faces
1225  for (unsigned i_face = 0; i_face < n_faces_attached_to_edge;
1226  i_face++)
1227  {
1228  // Is the node periodic in the face direction?
1229  is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
1230  faces_attached_to_edge[i_face]);
1231 
1232  // Check if the edge is periodic in the i_face-th face direction
1233  if (is_periodic)
1234  {
1235  // We're done!
1236  break;
1237  }
1238  } // for (unsigned
1239  // i_face=0;i_face<n_faces_attached_to_edge;i_face++)
1240 
1241  // Return the neighbour node pointer
1242  return neighbour_node_pt;
1243  } // if (neighbour_node_pt!=0)
1244  } // if (neigh_pt->object_pt()->nodes_built())
1245  } // if (neigh_pt!=0)
1246 
1247  // Keep searching, but only if there are further edge neighbours
1248  // Try next root edge neighbour
1249  i_root_edge_neighbour++;
1250 
1251  // Have we exhausted the search?
1252  if (i_root_edge_neighbour >= nroot_edge_neighbour)
1253  {
1254  // Stop searching
1255  keep_searching = false;
1256  }
1257  } // End of while keep searching over all true edge neighbours
1258  } // End of loop over edges
1259 
1260  // Node not found, return null
1261  return 0;
1262  }
1263 
1264 
1265  //==================================================================
1266  /// Build the element by doing the following:
1267  /// - Give it nodal positions (by establishing the pointers to its
1268  /// nodes)
1269  /// - In the process create new nodes where required (i.e. if they
1270  /// don't exist in father element or have already been created
1271  /// while building new neighbour elements). Node building
1272  /// involves the following steps:
1273  /// - Get nodal position from father element.
1274  /// - Establish the time-history of the newly created nodal point
1275  /// (its coordinates and the previous values) consistent with
1276  /// the father's history.
1277  /// - Determine the boundary conditions of the nodes (newly
1278  /// created nodes can only lie on the interior of any
1279  /// edges of the father element -- this makes it possible to
1280  /// to figure out what their bc should be...)
1281  /// - Add node to the mesh's stoarge scheme for the boundary nodes.
1282  /// - Add the new node to the mesh itself
1283  /// - Doc newly created nodes in file "new_nodes.dat" in the directory
1284  /// / of the DocInfo object (only if it's open!)
1285  /// - Finally, excute the element-specific further_build()
1286  /// (empty by default -- must be overloaded for specific elements).
1287  /// This deals with any build operations that are not included
1288  /// in the generic process outlined above. For instance, in
1289  /// Crouzeix Raviart elements we need to initialise the internal
1290  /// pressure values in manner consistent with the pressure
1291  /// distribution in the father element.
1292  //==================================================================
1293  void RefineableQElement<3>::build(Mesh*& mesh_pt,
1294  Vector<Node*>& new_node_pt,
1295  bool& was_already_built,
1296  std::ofstream& new_nodes_file)
1297  {
1298  using namespace OcTreeNames;
1299 
1300  // Number of dimensions
1301  unsigned n_dim = 3;
1302 
1303  // Get the number of 1d nodes
1304  unsigned n_p = nnode_1d();
1305 
1306  // Check whether static father_bound needs to be created
1307  if (Father_bound[n_p].nrow() == 0)
1308  {
1309  setup_father_bounds();
1310  }
1311 
1312  // Pointer to my father (in octree impersonation)
1313  OcTree* father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
1314 
1315  // What type of son am I? Ask my octree representation...
1316  int son_type = octree_pt()->son_type();
1317 
1318  // Has somebody build me already? (If any nodes have not been built)
1319  if (!nodes_built())
1320  {
1321 #ifdef PARANOID
1322  if (father_pt == 0)
1323  {
1324  std::string error_message =
1325  "Something fishy here: I have no father and yet \n";
1326  error_message += "I have no nodes. Who has created me then?!\n";
1327 
1328  throw OomphLibError(
1329  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1330  }
1331 #endif
1332 
1333  // Indicate status:
1334  was_already_built = false;
1335 
1336  // Return pointer to father element
1337  RefineableQElement<3>* father_el_pt =
1338  dynamic_cast<RefineableQElement<3>*>(father_pt->object_pt());
1339 
1340  // Timestepper should be the same for all nodes in father
1341  // element -- use it create timesteppers for new nodes
1342  TimeStepper* time_stepper_pt =
1343  father_el_pt->node_pt(0)->time_stepper_pt();
1344 
1345  // Number of history values (incl. present)
1346  unsigned ntstorage = time_stepper_pt->ntstorage();
1347 
1348  // Currently we can't handle the case of generalised coordinates
1349  // since we haven't established how they should be interpolated
1350  // Buffer this case:
1351  if (father_el_pt->node_pt(0)->nposition_type() != 1)
1352  {
1353  throw OomphLibError("Can't handle generalised nodal positions (yet).",
1354  OOMPH_CURRENT_FUNCTION,
1355  OOMPH_EXCEPTION_LOCATION);
1356  }
1357 
1358  Vector<int> s_lo(n_dim);
1359  Vector<int> s_hi(n_dim);
1360  Vector<double> s(n_dim);
1361  Vector<double> x(n_dim);
1362 
1363  // Setup vertex coordinates in father element:
1364  //--------------------------------------------
1365  // find the s_lo coordinates
1366  s_lo = octree_pt()->Direction_to_vector[son_type];
1367 
1368  // Just scale them, because the Direction_to_vector
1369  // doesn't really gives s_lo;
1370  for (unsigned i = 0; i < n_dim; i++)
1371  {
1372  s_lo[i] = (s_lo[i] + 1) / 2 - 1;
1373  }
1374 
1375  // setup s_hi (Actually s_hi[i]=s_lo[i]+1)
1376  for (unsigned i = 0; i < n_dim; i++)
1377  {
1378  s_hi[i] = s_lo[i] + 1;
1379  }
1380 
1381  // Pass macro element pointer on to sons and
1382  // set coordinates in macro element
1383  if (father_el_pt->macro_elem_pt() != 0)
1384  {
1385  set_macro_elem_pt(father_el_pt->macro_elem_pt());
1386  for (unsigned i = 0; i < n_dim; i++)
1387  {
1388  s_macro_ll(i) =
1389  father_el_pt->s_macro_ll(i) +
1390  0.5 * (s_lo[i] + 1.0) *
1391  (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1392  s_macro_ur(i) =
1393  father_el_pt->s_macro_ll(i) +
1394  0.5 * (s_hi[i] + 1.0) *
1395  (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1396  }
1397  }
1398 
1399 
1400  // If the father element hasn't been generated yet, we're stuck...
1401  if (father_el_pt->node_pt(0) == 0)
1402  {
1403  throw OomphLibError(
1404  "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1405  OOMPH_CURRENT_FUNCTION,
1406  OOMPH_EXCEPTION_LOCATION);
1407  }
1408  else
1409  {
1410  unsigned jnod = 0;
1411  Vector<double> s_fraction(n_dim);
1412 
1413  // Loop over nodes in element
1414  for (unsigned i0 = 0; i0 < n_p; i0++)
1415  {
1416  // Get the fractional position of the node in the direction of s[0]
1417  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
1418  // Local coordinate in father element
1419  s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
1420 
1421  for (unsigned i1 = 0; i1 < n_p; i1++)
1422  {
1423  // Get the fractional position of the node in the direction of s[1]
1424  s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
1425  // Local coordinate in father element
1426  s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
1427 
1428  for (unsigned i2 = 0; i2 < n_p; i2++)
1429  {
1430  // Get the fractional position of the node in the direction of
1431  // s[2]
1432  s_fraction[2] = local_one_d_fraction_of_node(i2, 2);
1433 
1434  // Local coordinate in father element
1435  s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
1436 
1437  // Local node number
1438  jnod = i0 + n_p * i1 + n_p * n_p * i2;
1439 
1440  // Initialise flag: So far, this node hasn't been built
1441  // or copied yet
1442  bool node_done = false;
1443 
1444  // Get the pointer to the node in the father; returns NULL
1445  // if there is not a node
1446  Node* created_node_pt =
1447  father_el_pt->get_node_at_local_coordinate(s);
1448 
1449  // Does this node already exist in father element?
1450  //------------------------------------------------
1451  bool node_exists_in_father = false;
1452  if (created_node_pt != 0)
1453  {
1454  // Remember this!
1455  node_exists_in_father = true;
1456 
1457  // Copy node across
1458  node_pt(jnod) = created_node_pt;
1459 
1460  // Make sure that we update the values at the node so that
1461  // they are consistent with the present representation.
1462  // This is only need for mixed interpolation where the value
1463  // at the father could now become active.
1464 
1465  // Loop over all history values
1466  for (unsigned t = 0; t < ntstorage; t++)
1467  {
1468  // Get values from father element
1469  // Note: get_interpolated_values() sets Vector size itself.
1470  Vector<double> prev_values;
1471  father_el_pt->get_interpolated_values(t, s, prev_values);
1472  // Find the minimum number of values
1473  //(either those stored at the node, or those returned by
1474  // the function)
1475  unsigned n_val_at_node = created_node_pt->nvalue();
1476  unsigned n_val_from_function = prev_values.size();
1477  // Use the ternary conditional operator here
1478  unsigned n_var = n_val_at_node < n_val_from_function ?
1479  n_val_at_node :
1480  n_val_from_function;
1481  // Assign the values that we can
1482  for (unsigned k = 0; k < n_var; k++)
1483  {
1484  created_node_pt->set_value(t, k, prev_values[k]);
1485  }
1486  }
1487 
1488  // Node has been created by copy
1489  node_done = true;
1490  }
1491  // Node does not exist in father element but might already
1492  //--------------------------------------------------------
1493  // have been created by neighbouring elements
1494  //-------------------------------------------
1495  else
1496  {
1497  // Boolean to check if the node is periodic
1498  bool is_periodic = false;
1499 
1500  // Was the node created by one of its neighbours
1501  // Whether or not the node lies on an edge can be determined
1502  // from the fractional position
1503  created_node_pt =
1504  node_created_by_neighbour(s_fraction, is_periodic);
1505 
1506  // If so, then copy the pointer across
1507  if (created_node_pt != 0)
1508  {
1509  // Now the node must be on a boundary, but we don't know which
1510  // one. The returned created_node_pt is actually the
1511  // neighbouring periodic node
1512  Node* neighbour_node_pt = created_node_pt;
1513 
1514  // Determine the edge on which the new node will live
1515  int father_bound = Father_bound[n_p](jnod, son_type);
1516 
1517  // Storage for the set of Mesh boundaries on which the
1518  // appropriate father edge lives.
1519  std::set<unsigned> boundaries;
1520 
1521  // Only get the boundaries if we are at the edge of
1522  // an element. Nodes in the centre of an element cannot be
1523  // on Mesh boundaries
1524  if (father_bound != Tree::OMEGA)
1525  {
1526  father_el_pt->get_boundaries(father_bound, boundaries);
1527  }
1528 
1529 #ifdef PARANOID
1530  // Case where a new node lives on more than one boundary
1531  // seems fishy enough to flag
1532  if (boundaries.size() > 2)
1533  {
1534  throw OomphLibError(
1535  "boundaries.size()>2 seems a bit strange..\n",
1536  OOMPH_CURRENT_FUNCTION,
1537  OOMPH_EXCEPTION_LOCATION);
1538  }
1539 #endif
1540 
1541  // If the node is periodic and is definitely a boundary node.
1542  // NOTE: the reason for this is that confusion can arise when
1543  // a node is created on an edge that joins a periodic face.
1544  if ((is_periodic) && (boundaries.size() > 0))
1545  {
1546  // Create node and set the pointer to it from the element
1547  created_node_pt =
1548  construct_boundary_node(jnod, time_stepper_pt);
1549 
1550  // Make the node periodic from the neighbour
1551  created_node_pt->make_periodic(neighbour_node_pt);
1552 
1553  // Add to vector of new nodes
1554  new_node_pt.push_back(created_node_pt);
1555 
1556  // Loop over # of history values
1557  for (unsigned t = 0; t < ntstorage; t++)
1558  {
1559  // Get position from father element -- this uses the macro
1560  // element representation if appropriate. If the node
1561  // turns out to be a hanging node later on, then
1562  // its position gets adjusted in line with its
1563  // hanging node interpolation.
1564  Vector<double> x_prev(n_dim);
1565  father_el_pt->get_x(t, s, x_prev);
1566  // Set previous positions of the new node
1567  for (unsigned i = 0; i < n_dim; i++)
1568  {
1569  created_node_pt->x(t, i) = x_prev[i];
1570  }
1571  }
1572 
1573  // Next, we Update the boundary lookup schemes
1574  // Loop over the boundaries stored in the set
1575  for (std::set<unsigned>::iterator it = boundaries.begin();
1576  it != boundaries.end();
1577  ++it)
1578  {
1579  // Add the node to the boundary
1580  mesh_pt->add_boundary_node(*it, created_node_pt);
1581 
1582  // If we have set an intrinsic coordinate on this
1583  // mesh boundary then it must also be interpolated on
1584  // the new node
1585  // Now interpolate the intrinsic boundary coordinate
1586  if (mesh_pt->boundary_coordinate_exists(*it) == true)
1587  {
1588  Vector<double> zeta(2, 0.0);
1589  father_el_pt->interpolated_zeta_on_face(
1590  *it, father_bound, s, zeta);
1591  created_node_pt->set_coordinates_on_boundary(*it, zeta);
1592  }
1593  }
1594 
1595  // Make sure that we add the node to the mesh
1596  mesh_pt->add_node_pt(created_node_pt);
1597  } // End of periodic case
1598  // Otherwise the node is not periodic, so just set the
1599  // pointer to the neighbours node
1600  else
1601  {
1602  node_pt(jnod) = created_node_pt;
1603  }
1604  node_done = true;
1605  }
1606  // Node does not exist in neighbour element but might already
1607  //-----------------------------------------------------------
1608  // have been created by a son of a neighbouring element
1609  //-----------------------------------------------------
1610  else
1611  {
1612  // Was the node created by one of its neighbours' sons
1613  // Whether or not the node lies on an edge can be calculated
1614  // by from the fractional position
1615  bool is_periodic = false;
1616  created_node_pt =
1617  node_created_by_son_of_neighbour(s_fraction, is_periodic);
1618 
1619  // If the node was so created, assign the pointers
1620  if (created_node_pt != 0)
1621  {
1622  // If the node is periodic
1623  if (is_periodic)
1624  {
1625  // Now the node must be on a boundary, but we don't know
1626  // which one The returned created_node_pt is actually the
1627  // neighbouring periodic node
1628  Node* neighbour_node_pt = created_node_pt;
1629 
1630  // Determine the edge on which the new node will live
1631  int father_bound = Father_bound[n_p](jnod, son_type);
1632 
1633  // Storage for the set of Mesh boundaries on which the
1634  // appropriate father edge lives.
1635  std::set<unsigned> boundaries;
1636 
1637  // Only get the boundaries if we are at the edge of
1638  // an element. Nodes in the centre of an element cannot be
1639  // on Mesh boundaries
1640  if (father_bound != Tree::OMEGA)
1641  {
1642  father_el_pt->get_boundaries(father_bound, boundaries);
1643  }
1644 
1645 #ifdef PARANOID
1646  // Case where a new node lives on more than one boundary
1647  // seems fishy enough to flag
1648  if (boundaries.size() > 2)
1649  {
1650  throw OomphLibError(
1651  "boundaries.size()>2 seems a bit strange..\n",
1652  OOMPH_CURRENT_FUNCTION,
1653  OOMPH_EXCEPTION_LOCATION);
1654  }
1655 
1656  // Case when there are no boundaries, we are in big
1657  // trouble
1658  if (boundaries.size() == 0)
1659  {
1660  std::ostringstream error_stream;
1661  error_stream
1662  << "Periodic node is not on a boundary...\n"
1663  << "Coordinates: " << created_node_pt->x(0) << " "
1664  << created_node_pt->x(1) << "\n";
1665  throw OomphLibError(error_stream.str(),
1666  OOMPH_CURRENT_FUNCTION,
1667  OOMPH_EXCEPTION_LOCATION);
1668  }
1669 #endif
1670 
1671  // Create node and set the pointer to it from the element
1672  created_node_pt =
1673  construct_boundary_node(jnod, time_stepper_pt);
1674 
1675  // Make the node periodic from the neighbour
1676  created_node_pt->make_periodic(neighbour_node_pt);
1677 
1678  // Add to vector of new nodes
1679  new_node_pt.push_back(created_node_pt);
1680 
1681  // Loop over # of history values
1682  for (unsigned t = 0; t < ntstorage; t++)
1683  {
1684  // Get position from father element -- this uses the
1685  // macro element representation if appropriate. If the
1686  // node turns out to be a hanging node later on, then
1687  // its position gets adjusted in line with its
1688  // hanging node interpolation.
1689  Vector<double> x_prev(n_dim, 0.0);
1690  father_el_pt->get_x(t, s, x_prev);
1691  // Set previous positions of the new node
1692  for (unsigned i = 0; i < n_dim; i++)
1693  {
1694  created_node_pt->x(t, i) = x_prev[i];
1695  }
1696  }
1697 
1698  // Next, we Update the boundary lookup schemes
1699  // Loop over the boundaries stored in the set
1700  for (std::set<unsigned>::iterator it = boundaries.begin();
1701  it != boundaries.end();
1702  ++it)
1703  {
1704  // Add the node to the boundary
1705  mesh_pt->add_boundary_node(*it, created_node_pt);
1706 
1707  // If we have set an intrinsic coordinate on this
1708  // mesh boundary then it must also be interpolated on
1709  // the new node
1710  // Now interpolate the intrinsic boundary coordinate
1711  if (mesh_pt->boundary_coordinate_exists(*it) == true)
1712  {
1713  Vector<double> zeta(2, 0.0);
1714  father_el_pt->interpolated_zeta_on_face(
1715  *it, father_bound, s, zeta);
1716  created_node_pt->set_coordinates_on_boundary(*it,
1717  zeta);
1718  }
1719  }
1720 
1721  // Make sure that we add the node to the mesh
1722  mesh_pt->add_node_pt(created_node_pt);
1723  } // End of periodic case
1724  // Otherwise the node is not periodic, so just set the
1725  // pointer to the neighbours node
1726  else
1727  {
1728  node_pt(jnod) = created_node_pt;
1729  }
1730  // Node has been created
1731  node_done = true;
1732  } // Node does not exist in son of neighbouring element
1733  } // Node does not exist in neighbouring element
1734  } // Node does not exist in father element
1735 
1736  // Node already exists in father: No need to do anything else!
1737  // otherwise deal with boundary information etc.
1738  if (!node_exists_in_father)
1739  {
1740  // Check boundary status
1741  //----------------------
1742 
1743  // Firstly, we need to determine whether or not a node lies
1744  // on the boundary before building it, because
1745  // we actually assign a different type of node on boundaries.
1746 
1747  // If the new node lives on a face that is
1748  // shared with a face of its father element,
1749  // it needs to inherit the bounday conditions
1750  // from the father face
1751  int father_bound = Father_bound[n_p](jnod, son_type);
1752 
1753  // Storage for the set of Mesh boundaries on which the
1754  // appropriate father face lives.
1755  // [New nodes should always be mid-edge nodes in father
1756  // and therefore only live on one boundary but just to
1757  // play it safe...]
1758  std::set<unsigned> boundaries;
1759 
1760  // Only get the boundaries if we are at the edge of
1761  // an element. Nodes in the centre of an element cannot be
1762  // on Mesh boundaries
1763  if (father_bound != Tree::OMEGA)
1764  {
1765  father_el_pt->get_boundaries(father_bound, boundaries);
1766  }
1767 
1768 #ifdef PARANOID
1769  // Case where a new node lives on more than two boundaries
1770  // seems fishy enough to flag
1771  if (boundaries.size() > 2)
1772  {
1773  throw OomphLibError(
1774  "boundaries.size()>2 seems a bit strange..\n",
1775  OOMPH_CURRENT_FUNCTION,
1776  OOMPH_EXCEPTION_LOCATION);
1777  }
1778 #endif
1779 
1780  // If the node lives on a mesh boundary,
1781  // then we need to create a boundary node
1782  if (boundaries.size() > 0)
1783  {
1784  // Do we need a new node?
1785  if (!node_done)
1786  {
1787  // Create node and set the internal pointer
1788  created_node_pt =
1789  construct_boundary_node(jnod, time_stepper_pt);
1790  // Add to vector of new nodes
1791  new_node_pt.push_back(created_node_pt);
1792  }
1793 
1794  // Now we need to work out whether to pin the values at
1795  // the new node based on the boundary conditions applied at
1796  // its Mesh boundary
1797 
1798  // Get the boundary conditions from the father.
1799  // Note: We can only deal with the values that are
1800  // continuously interpolated in the bulk element.
1801  unsigned n_cont = ncont_interpolated_values();
1802  Vector<int> bound_cons(n_cont);
1803  father_el_pt->get_bcs(father_bound, bound_cons);
1804 
1805  // Loop over the continuously interpolated values and pin,
1806  // if necessary
1807  for (unsigned k = 0; k < n_cont; k++)
1808  {
1809  if (bound_cons[k])
1810  {
1811  created_node_pt->pin(k);
1812  }
1813  }
1814 
1815  // Solid node? If so, deal with the positional boundary
1816  // conditions:
1817  SolidNode* solid_node_pt =
1818  dynamic_cast<SolidNode*>(created_node_pt);
1819  if (solid_node_pt != 0)
1820  {
1821  // Get the positional boundary conditions from the father:
1822  unsigned n_dim = created_node_pt->ndim();
1823  Vector<int> solid_bound_cons(n_dim);
1824  RefineableSolidQElement<3>* father_solid_el_pt =
1825  dynamic_cast<RefineableSolidQElement<3>*>(father_el_pt);
1826 #ifdef PARANOID
1827  if (father_solid_el_pt == 0)
1828  {
1829  std::string error_message = "We have a SolidNode outside "
1830  "a refineable SolidElement\n";
1831  error_message +=
1832  "during mesh refinement -- this doesn't make sense";
1833 
1834  throw OomphLibError(error_message,
1835  OOMPH_CURRENT_FUNCTION,
1836  OOMPH_EXCEPTION_LOCATION);
1837  }
1838 #endif
1839  father_solid_el_pt->get_solid_bcs(father_bound,
1840  solid_bound_cons);
1841 
1842  // Loop over the positions and pin, if necessary
1843  for (unsigned k = 0; k < n_dim; k++)
1844  {
1845  if (solid_bound_cons[k])
1846  {
1847  solid_node_pt->pin_position(k);
1848  }
1849  }
1850  } // End of if solid_node_pt
1851 
1852  // Next update the boundary look-up schemes
1853 
1854  // Loop over the boundaries in the set
1855  for (std::set<unsigned>::iterator it = boundaries.begin();
1856  it != boundaries.end();
1857  ++it)
1858  {
1859  // Add the node to the bounadry
1860  mesh_pt->add_boundary_node(*it, created_node_pt);
1861 
1862  // If we have set an intrinsic coordinate on this
1863  // mesh boundary then it must also be interpolated on
1864  // the new node
1865  // Now interpolate the intrinsic boundary coordinate
1866  if (mesh_pt->boundary_coordinate_exists(*it) == true)
1867  {
1868  // Usually there will be two coordinates
1869  Vector<double> zeta(2);
1870  father_el_pt->interpolated_zeta_on_face(
1871  *it, father_bound, s, zeta);
1872 
1873  created_node_pt->set_coordinates_on_boundary(*it, zeta);
1874  }
1875  }
1876  }
1877  // Otherwise the node is not on a Mesh boundary and
1878  // we create a normal "bulk" node
1879  else
1880  {
1881  // Do we need a new node?
1882  if (!node_done)
1883  {
1884  // Create node and set the pointer to it from the element
1885  created_node_pt = construct_node(jnod, time_stepper_pt);
1886  // Add to vector of new nodes
1887  new_node_pt.push_back(created_node_pt);
1888  }
1889  }
1890 
1891  // In the first instance use macro element or FE representation
1892  // to create past and present nodal positions.
1893  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
1894  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
1895  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
1896  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
1897  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
1898  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
1899  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
1900 
1901 
1902  // Have we created a new node?
1903  if (!node_done)
1904  {
1905  // Loop over # of history values
1906  for (unsigned t = 0; t < ntstorage; t++)
1907  {
1908  // Get position from father element -- this uses the macro
1909  // element representation if appropriate. If the node
1910  // turns out to be a hanging node later on, then
1911  // its position gets adjusted in line with its
1912  // hanging node interpolation.
1913  Vector<double> x_prev(n_dim);
1914  father_el_pt->get_x(t, s, x_prev);
1915 
1916  // Set previous positions of the new node
1917  for (unsigned i = 0; i < n_dim; i++)
1918  {
1919  created_node_pt->x(t, i) = x_prev[i];
1920  }
1921  }
1922 
1923  // Now set the values
1924  // Loop over all history values
1925  for (unsigned t = 0; t < ntstorage; t++)
1926  {
1927  // Get values from father element
1928  // Note: get_interpolated_values() sets Vector size itself.
1929  Vector<double> prev_values;
1930  father_el_pt->get_interpolated_values(t, s, prev_values);
1931 
1932  // Initialise the values at the new node
1933  unsigned n_value = created_node_pt->nvalue();
1934  for (unsigned k = 0; k < n_value; k++)
1935  {
1936  created_node_pt->set_value(t, k, prev_values[k]);
1937  }
1938  }
1939 
1940  // Add new node to mesh
1941  mesh_pt->add_node_pt(created_node_pt);
1942 
1943  } // End of whether the node has been created by us or not
1944 
1945  } // End of if node already existed in father in which case
1946  // everything above gets bypassed.
1947 
1948  // Check if the element is an algebraic element
1949  AlgebraicElementBase* alg_el_pt =
1950  dynamic_cast<AlgebraicElementBase*>(this);
1951 
1952  // If the element is an algebraic element, setup
1953  // node position (past and present) from algebraic update
1954  // function.
1955  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
1956  // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
1957  // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
1958  // INFO FOR *ALL* ROOT ELEMENTS!
1959  if (alg_el_pt != 0)
1960  {
1961  // Build algebraic node update info for new node
1962  // This sets up the node update data for all node update
1963  // functions that are shared by all nodes in the father
1964  // element
1965  alg_el_pt->setup_algebraic_node_update(
1966  node_pt(jnod), s, father_el_pt);
1967  }
1968 
1969  // We have built the node and we are documenting
1970  if ((!node_done) && (new_nodes_file.is_open()))
1971  {
1972  new_nodes_file << node_pt(jnod)->x(0) << " "
1973  << node_pt(jnod)->x(1) << " "
1974  << node_pt(jnod)->x(2) << std::endl;
1975  }
1976 
1977  } // End of Z loop over nodes in element
1978 
1979  } // End of vertical loop over nodes in element
1980 
1981  } // End of horizontal loop over nodes in element
1982 
1983  // If the element is a MacroElementNodeUpdateElement, set
1984  // the update parameters for the current element's nodes --
1985  // all this needs is the vector of (pointers to the)
1986  // geometric objects that affect the MacroElement-based
1987  // node update -- this is the same as that in the father element
1988  MacroElementNodeUpdateElementBase* father_m_el_pt =
1989  dynamic_cast<MacroElementNodeUpdateElementBase*>(father_el_pt);
1990  if (father_m_el_pt != 0)
1991  {
1992  // Get vector of geometric objects from father (construct vector
1993  // via copy operation)
1994  Vector<GeomObject*> geom_object_pt(father_m_el_pt->geom_object_pt());
1995 
1996  // Cast current element to MacroElementNodeUpdateElement:
1997  MacroElementNodeUpdateElementBase* m_el_pt =
1998  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
1999 
2000 #ifdef PARANOID
2001  if (m_el_pt == 0)
2002  {
2003  std::string error_message =
2004  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
2005  error_message +=
2006  "Strange -- if the father is a MacroElementNodeUpdateElement\n";
2007  error_message += "the son should be too....\n";
2008 
2009  throw OomphLibError(
2010  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2011  }
2012 #endif
2013  // Build update info by passing vector of geometric objects:
2014  // This sets the current element to be the update element
2015  // for all of the element's nodes -- this is reversed
2016  // if the element is ever un-refined in the father element's
2017  // rebuild_from_sons() function which overwrites this
2018  // assignment to avoid nasty segmentation faults that occur
2019  // when a node tries to update itself via an element that no
2020  // longer exists...
2021  m_el_pt->set_node_update_info(geom_object_pt);
2022  }
2023 
2024 #ifdef OOMPH_HAS_MPI
2025  // Is the new element a halo element?
2026  if (tree_pt()->father_pt()->object_pt()->is_halo())
2027  {
2028  Non_halo_proc_ID =
2029  tree_pt()->father_pt()->object_pt()->non_halo_proc_ID();
2030  }
2031 #endif
2032 
2033  // Is it an ElementWithMovingNodes?
2034  ElementWithMovingNodes* aux_el_pt =
2035  dynamic_cast<ElementWithMovingNodes*>(this);
2036 
2037  // Pass down the information re the method for the evaluation
2038  // of the shape derivatives
2039  if (aux_el_pt != 0)
2040  {
2041  ElementWithMovingNodes* aux_father_el_pt =
2042  dynamic_cast<ElementWithMovingNodes*>(father_el_pt);
2043 
2044 #ifdef PARANOID
2045  if (aux_father_el_pt == 0)
2046  {
2047  std::string error_message =
2048  "Failed to cast to ElementWithMovingNodes*\n";
2049  error_message +=
2050  "Strange -- if the son is a ElementWithMovingNodes\n";
2051  error_message += "the father should be too....\n";
2052 
2053  throw OomphLibError(
2054  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2055  }
2056 #endif
2057 
2058  // If evaluating the residuals by finite differences in the father
2059  // continue to do so in the child
2060  if (aux_father_el_pt
2061  ->are_dresidual_dnodal_coordinates_always_evaluated_by_fd())
2062  {
2063  aux_el_pt
2064  ->enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
2065  }
2066 
2067 
2068  aux_el_pt->method_for_shape_derivs() =
2069  aux_father_el_pt->method_for_shape_derivs();
2070 
2071  // If bypassing the evaluation of fill_in_jacobian_from_geometric_data
2072  // continue to do so
2073  if (aux_father_el_pt
2074  ->is_fill_in_jacobian_from_geometric_data_bypassed())
2075  {
2076  aux_el_pt->enable_bypass_fill_in_jacobian_from_geometric_data();
2077  }
2078  }
2079 
2080  // Now do further build (if any)
2081  further_build();
2082 
2083  } // Sanity check: Father element has been generated
2084 
2085  } // End for element has not been built yet
2086  else
2087  {
2088  was_already_built = true;
2089  }
2090  }
2091 
2092  //====================================================================
2093  /// Set up all hanging nodes.
2094  ///
2095  /// (Wrapper to avoid specification of the unwanted output file).
2096  //====================================================================
2097  void RefineableQElement<3>::setup_hanging_nodes(
2098  Vector<std::ofstream*>& output_stream)
2099  {
2100 #ifdef PARANOID
2101  if (output_stream.size() != 6)
2102  {
2103  throw OomphLibError("There must be six output streams",
2104  OOMPH_CURRENT_FUNCTION,
2105  OOMPH_EXCEPTION_LOCATION);
2106  }
2107 #endif
2108 
2109  using namespace OcTreeNames;
2110 
2111  // Setup hanging nodes on each edge of the element
2112  oc_hang_helper(-1, U, *(output_stream[0]));
2113  oc_hang_helper(-1, D, *(output_stream[1]));
2114  oc_hang_helper(-1, L, *(output_stream[2]));
2115  oc_hang_helper(-1, R, *(output_stream[3]));
2116  oc_hang_helper(-1, B, *(output_stream[4]));
2117  oc_hang_helper(-1, F, *(output_stream[5]));
2118  }
2119 
2120  //================================================================
2121  /// Internal function that sets up the hanging node scheme for a
2122  /// particular value
2123  //=================================================================
2124  void RefineableQElement<3>::setup_hang_for_value(const int& value_id)
2125  {
2126  using namespace OcTreeNames;
2127 
2128  std::ofstream dummy_hangfile;
2129  // Setup hanging nodes on each edge of the element
2130  oc_hang_helper(value_id, U, dummy_hangfile);
2131  oc_hang_helper(value_id, D, dummy_hangfile);
2132  oc_hang_helper(value_id, R, dummy_hangfile);
2133  oc_hang_helper(value_id, L, dummy_hangfile);
2134  oc_hang_helper(value_id, B, dummy_hangfile);
2135  oc_hang_helper(value_id, F, dummy_hangfile);
2136  }
2137 
2138  //=================================================================
2139  /// Internal function to set up the hanging nodes on a particular
2140  /// face of the element
2141  //=================================================================
2142  void RefineableQElement<3>::oc_hang_helper(const int& value_id,
2143  const int& my_face,
2144  std::ofstream& output_hangfile)
2145  {
2146  using namespace OcTreeNames;
2147 
2148  // Number of dimensions
2149  unsigned n_dim = 3;
2150 
2151  Vector<unsigned> translate_s(n_dim);
2152  Vector<double> s_lo_neigh(n_dim);
2153  Vector<double> s_hi_neigh(n_dim);
2154  int neigh_face;
2155  int diff_level;
2156  bool in_neighbouring_tree;
2157 
2158  // Find pointer to neighbour in this direction
2159  OcTree* neigh_pt;
2160  neigh_pt = octree_pt()->gteq_face_neighbour(my_face,
2161  translate_s,
2162  s_lo_neigh,
2163  s_hi_neigh,
2164  neigh_face,
2165  diff_level,
2166  in_neighbouring_tree);
2167 
2168  // Neighbour exists
2169  if (neigh_pt != 0)
2170  {
2171  // Different sized element?
2172  if (diff_level != 0)
2173  {
2174  // Test for the periodic node case
2175  // Are we crossing a periodic boundary
2176  bool is_periodic = false;
2177  if (in_neighbouring_tree)
2178  {
2179  is_periodic = tree_pt()->root_pt()->is_neighbour_periodic(my_face);
2180  }
2181 
2182  // If it is periodic we actually need to get the node in
2183  // the neighbour of the neighbour (which will be a parent of
2184  // the present element) so that the "fixed" coordinate is
2185  // correctly calculated.
2186  // The idea is to replace the neigh_pt and associated data
2187  // with those of the neighbour of the neighbour
2188  if (is_periodic)
2189  {
2190  // Required data for the neighbour finding routine
2191  Vector<unsigned> translate_s_in_neigh(n_dim, 0.0);
2192  Vector<double> s_lo_neigh_of_neigh(n_dim, 0.0);
2193  Vector<double> s_hi_neigh_of_neigh(n_dim, 0.0);
2194  int neigh_face_of_neigh;
2195  int diff_level_of_neigh;
2196  bool in_neighbouring_tree_of_neigh;
2197 
2198  // Find pointer to neighbour of the neighbour on the edge
2199  // that we are currently considering
2200  OcTree* neigh_of_neigh_pt;
2201  neigh_of_neigh_pt =
2202  neigh_pt->gteq_face_neighbour(neigh_face,
2203  translate_s_in_neigh,
2204  s_lo_neigh_of_neigh,
2205  s_hi_neigh_of_neigh,
2206  neigh_face_of_neigh,
2207  diff_level_of_neigh,
2208  in_neighbouring_tree_of_neigh);
2209 
2210  // Set the value of the NEW neighbour and edge
2211  neigh_pt = neigh_of_neigh_pt;
2212  neigh_face = neigh_face_of_neigh;
2213 
2214  // Set the values of the translation schemes
2215  // Need to find the values of s_lo and s_hi
2216  // in the neighbour of the neighbour
2217 
2218  // Get the minimum and maximum values of the coordinate
2219  // in the neighbour (don't like this, but I think it's
2220  // necessary) Note that these values are hardcoded
2221  // in the quadtrees at some point!!
2222  double s_min = neigh_pt->object_pt()->s_min();
2223  double s_max = neigh_pt->object_pt()->s_max();
2224  Vector<double> s_lo_frac(n_dim), s_hi_frac(n_dim);
2225  // Work out the fractional position of the low and high points
2226  // of the original element
2227  for (unsigned i = 0; i < n_dim; i++)
2228  {
2229  s_lo_frac[i] = (s_lo_neigh[i] - s_min) / (s_max - s_min);
2230  s_hi_frac[i] = (s_hi_neigh[i] - s_min) / (s_max - s_min);
2231  }
2232 
2233  // We should now be able to construct the low and high points in
2234  // the neighbour of the neighbour
2235  for (unsigned i = 0; i < n_dim; i++)
2236  {
2237  s_lo_neigh[i] = s_lo_neigh_of_neigh[i] +
2238  s_lo_frac[translate_s_in_neigh[i]] *
2239  (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
2240  s_hi_neigh[i] = s_lo_neigh_of_neigh[i] +
2241  s_hi_frac[translate_s_in_neigh[i]] *
2242  (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
2243  }
2244 
2245  // Finally we must sort out the translation scheme
2246  Vector<unsigned> temp_translate(n_dim, 0.0);
2247  for (unsigned i = 0; i < n_dim; i++)
2248  {
2249  temp_translate[i] = translate_s_in_neigh[translate_s[i]];
2250  }
2251  for (unsigned i = 0; i < n_dim; i++)
2252  {
2253  translate_s[i] = temp_translate[i];
2254  }
2255  } // End of special treatment for periodic hanging nodes
2256 
2257  // Number of nodes in one dimension
2258  unsigned n_p = ninterpolating_node_1d(value_id);
2259  // Storage for the local nodes along the face of the element
2260  Node* local_node_pt = 0;
2261 
2262  // Loop over nodes along the face
2263  for (unsigned i0 = 0; i0 < n_p; i0++)
2264  {
2265  for (unsigned i1 = 0; i1 < n_p; i1++)
2266  {
2267  // Storage for the fractional position of the node in the element
2268  Vector<double> s_fraction(n_dim);
2269 
2270  // Local node number
2271  switch (my_face)
2272  {
2273  case U:
2274  s_fraction[0] =
2275  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2276  s_fraction[1] = 1.0;
2277  s_fraction[2] =
2278  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2279  local_node_pt = interpolating_node_pt(
2280  i0 + (n_p - 1) * n_p + n_p * n_p * i1, value_id);
2281  break;
2282 
2283  case D:
2284  s_fraction[0] =
2285  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2286  s_fraction[1] = 0.0;
2287  s_fraction[2] =
2288  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2289  local_node_pt =
2290  interpolating_node_pt(i0 + n_p * n_p * i1, value_id);
2291  break;
2292 
2293  case R:
2294  s_fraction[0] = 1.0;
2295  s_fraction[1] =
2296  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
2297  s_fraction[2] =
2298  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2299  local_node_pt = interpolating_node_pt(
2300  n_p - 1 + i0 * n_p + i1 * n_p * n_p, value_id);
2301  break;
2302 
2303  case L:
2304  s_fraction[0] = 0.0;
2305  s_fraction[1] =
2306  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
2307  s_fraction[2] =
2308  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2309  local_node_pt =
2310  interpolating_node_pt(n_p * i0 + i1 * n_p * n_p, value_id);
2311  break;
2312 
2313  case B:
2314  s_fraction[0] =
2315  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2316  s_fraction[1] =
2317  local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
2318  s_fraction[2] = 0.0;
2319  local_node_pt = interpolating_node_pt(i0 + i1 * n_p, value_id);
2320  break;
2321 
2322  case F:
2323  s_fraction[0] =
2324  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2325  s_fraction[1] =
2326  local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
2327  s_fraction[2] = 1.0;
2328  local_node_pt = interpolating_node_pt(
2329  i0 + i1 * n_p + (n_p - 1) * n_p * n_p, value_id);
2330  break;
2331 
2332  default:
2333  throw OomphLibError("my_face not U, D, L, R, B, F\n",
2334  OOMPH_CURRENT_FUNCTION,
2335  OOMPH_EXCEPTION_LOCATION);
2336  }
2337 
2338  // Set up the coordinate in the neighbour
2339  Vector<double> s_in_neighb(n_dim);
2340  for (unsigned i = 0; i < n_dim; i++)
2341  {
2342  s_in_neighb[i] =
2343  s_lo_neigh[i] +
2344  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
2345  }
2346 
2347  // Find the Node in the neighbouring element
2348  Node* neighbouring_node_pt =
2349  neigh_pt->object_pt()->get_interpolating_node_at_local_coordinate(
2350  s_in_neighb, value_id);
2351 
2352  // If the neighbour does not have a node at this point
2353  if (0 == neighbouring_node_pt)
2354  {
2355  // Do we need to make a hanging node, we assume that we don't
2356  // initially
2357  bool make_hanging_node = false;
2358 
2359  // If the node is not hanging geometrically, then we must make it
2360  // hang
2361  if (!local_node_pt->is_hanging())
2362  {
2363  make_hanging_node = true;
2364  }
2365  // Otherwise is could be hanging geometrically, but still require
2366  // a different hanging scheme for this data value
2367  else
2368  {
2369  if ((value_id != -1) && (local_node_pt->hanging_pt(value_id) ==
2370  local_node_pt->hanging_pt()))
2371  {
2372  make_hanging_node = true;
2373  }
2374  }
2375 
2376  if (make_hanging_node == true)
2377  {
2378  // Cache refineable element used here
2379  RefineableElement* const obj_pt = neigh_pt->object_pt();
2380 
2381  // Get shape functions in neighbour element
2382  Shape psi(obj_pt->ninterpolating_node(value_id));
2383  obj_pt->interpolating_basis(s_in_neighb, psi, value_id);
2384 
2385  // Allocate the storage for the Hang pointer
2386  // We know that it will hold n_p*n_p nodes
2387  HangInfo* hang_pt = new HangInfo(n_p * n_p);
2388 
2389  // Loop over nodes on edge in neighbour and mark them as nodes
2390  // that this node depends on
2391  unsigned n_neighbour;
2392 
2393  // Number of nodes along edge in neighbour element
2394  for (unsigned ii0 = 0; ii0 < n_p; ii0++)
2395  {
2396  for (unsigned ii1 = 0; ii1 < n_p; ii1++)
2397  {
2398  switch (neigh_face)
2399  {
2400  case U:
2401  n_neighbour = ii0 + n_p * (n_p - 1) + ii1 * n_p * n_p;
2402  break;
2403 
2404  case D:
2405  n_neighbour = ii0 + ii1 * n_p * n_p;
2406  break;
2407 
2408  case L:
2409  n_neighbour = ii0 * n_p + ii1 * n_p * n_p;
2410  break;
2411 
2412  case R:
2413  n_neighbour = (n_p - 1) + ii0 * n_p + ii1 * n_p * n_p;
2414  break;
2415 
2416  case B:
2417  n_neighbour = ii0 + ii1 * n_p;
2418  break;
2419 
2420  case F:
2421  n_neighbour = ii0 + ii1 * n_p + n_p * n_p * (n_p - 1);
2422  break;
2423  default:
2424  throw OomphLibError("neigh_face not U, L, R, B, F\n",
2425  OOMPH_CURRENT_FUNCTION,
2426  OOMPH_EXCEPTION_LOCATION);
2427  }
2428 
2429  // Push back neighbouring node and weight into
2430  // Vector of (pointers to)
2431  // master nodes and weights
2432  // The weight is merely the value of the shape function
2433  // corresponding to the node in the neighbour
2434  hang_pt->set_master_node_pt(
2435  ii0 * n_p + ii1,
2436  obj_pt->interpolating_node_pt(n_neighbour, value_id),
2437  psi[n_neighbour]);
2438  }
2439  }
2440  // Now set the hanging data for the position
2441  // This also constrains the data values associated with the
2442  // value id
2443  local_node_pt->set_hanging_pt(hang_pt, value_id);
2444  }
2445 
2446  if (output_hangfile.is_open())
2447  {
2448  // output_hangfile
2449  output_hangfile << local_node_pt->x(0) << " "
2450  << local_node_pt->x(1) << " "
2451  << local_node_pt->x(2) << std::endl;
2452  }
2453  }
2454  else
2455  {
2456 #ifdef PARANOID
2457  if (local_node_pt != neighbouring_node_pt)
2458  {
2459  std::ofstream reportage("dodgy.dat", std::ios_base::app);
2460  reportage << local_node_pt->x(0) << " " << local_node_pt->x(1)
2461  << " " << local_node_pt->x(2) << std::endl;
2462  reportage.close();
2463 
2464  std::ostringstream warning_stream;
2465  warning_stream
2466  << "SANITY CHECK in oc_hang_helper \n"
2467  << "Current node " << local_node_pt << " at "
2468  << "(" << local_node_pt->x(0) << ", " << local_node_pt->x(1)
2469  << ", " << local_node_pt->x(2) << ")" << std::endl
2470  << " is not hanging and has " << std::endl
2471  << "Neighbour's node " << neighbouring_node_pt << " at "
2472  << "(" << neighbouring_node_pt->x(0) << ", "
2473  << neighbouring_node_pt->x(1) << ", "
2474  << neighbouring_node_pt->x(2) << ")" << std::endl
2475  << "even though the two should be "
2476  << "identical" << std::endl;
2477  OomphLibWarning(warning_stream.str(),
2478  OOMPH_CURRENT_FUNCTION,
2479  OOMPH_EXCEPTION_LOCATION);
2480  }
2481 #endif
2482  }
2483 
2484  // If we are doing the position then
2485  if (value_id == -1)
2486  {
2487  // Get the nodal position from neighbour element
2488  Vector<double> x_in_neighb(n_dim);
2489  neigh_pt->object_pt()->interpolated_x(s_in_neighb, x_in_neighb);
2490 
2491  // Fine adjust the coordinates (macro map will pick up boundary
2492  // accurately but will lead to different element edges)
2493  local_node_pt->x(0) = x_in_neighb[0];
2494  local_node_pt->x(1) = x_in_neighb[1];
2495  local_node_pt->x(2) = x_in_neighb[2];
2496  }
2497  }
2498  }
2499  }
2500  }
2501  }
2502 
2503 
2504  //=================================================================
2505  /// Check inter-element continuity of
2506  /// - nodal positions
2507  /// - (nodally) interpolated function values
2508  //====================================================================
2509  void RefineableQElement<3>::check_integrity(double& max_error)
2510  {
2511  using namespace OcTreeNames;
2512 
2513  // std::ofstream error_file("errors.dat");
2514 
2515  // Number of nodes along edge
2516  unsigned n_p = nnode_1d();
2517 
2518  // Number of timesteps (incl. present) for which continuity is
2519  // to be checked.
2520  unsigned n_time = 1;
2521 
2522  // Initialise errors
2523  max_error = 0.0;
2524  Vector<double> max_error_x(3, 0.0);
2525  double max_error_val = 0.0;
2526 
2527  // Set the faces
2528  Vector<int> faces(6);
2529  faces[0] = D;
2530  faces[1] = U;
2531  faces[2] = L;
2532  faces[3] = R;
2533  faces[4] = B;
2534  faces[5] = F;
2535 
2536  // Loop over the edges
2537  for (unsigned face_counter = 0; face_counter < 6; face_counter++)
2538  {
2539  Vector<double> s(3), s_lo_neigh(3), s_hi_neigh(3);
2540  Vector<double> s_fraction(3);
2541  Vector<unsigned> translate_s(3);
2542  int neigh_face, diff_level;
2543  int my_face;
2544  bool in_neighbouring_tree;
2545 
2546  my_face = faces[face_counter];
2547 
2548  // Find pointer to neighbour in this direction
2549  OcTree* neigh_pt;
2550  neigh_pt = octree_pt()->gteq_face_neighbour(my_face,
2551  translate_s,
2552  s_lo_neigh,
2553  s_hi_neigh,
2554  neigh_face,
2555  diff_level,
2556  in_neighbouring_tree);
2557 
2558  // Neighbour exists and has existing nodes
2559  if ((neigh_pt != 0) && (neigh_pt->object_pt()->nodes_built()))
2560  {
2561  // if (diff_level!=0)
2562  {
2563  // Need to exclude periodic nodes from this check
2564  // There are only periodic nodes if we are in a neighbouring tree
2565  bool is_periodic = false;
2566  if (in_neighbouring_tree)
2567  {
2568  // Is it periodic
2569  is_periodic =
2570  tree_pt()->root_pt()->is_neighbour_periodic(faces[face_counter]);
2571  }
2572 
2573  // Loop over nodes along the edge
2574  for (unsigned i0 = 0; i0 < n_p; i0++)
2575  {
2576  for (unsigned i1 = 0; i1 < n_p; i1++)
2577  {
2578  // Local node number
2579  unsigned n = 0;
2580  switch (face_counter)
2581  {
2582  case 0:
2583  // Fractions
2584  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2585  s_fraction[1] = 0.0;
2586  s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2587  // Set local node number
2588  n = i0 + i1 * n_p * n_p;
2589  break;
2590 
2591  case 1:
2592  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2593  s_fraction[1] = 1.0;
2594  s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2595  // Set local node number
2596  n = i0 + n_p * (n_p - 1) + i1 * n_p * n_p;
2597  break;
2598 
2599  case 2:
2600  s_fraction[0] = 0.0;
2601  s_fraction[1] = local_one_d_fraction_of_node(i0, 1);
2602  s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2603  // Set local node number
2604  n = n_p * i0 + i1 * n_p * n_p;
2605  break;
2606 
2607  case 3:
2608  s_fraction[0] = 1.0;
2609  s_fraction[1] = local_one_d_fraction_of_node(i0, 1);
2610  s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2611  // Set local node number
2612  n = n_p - 1 + n_p * i0 + n_p * n_p * i1;
2613  break;
2614 
2615  case 4:
2616  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2617  s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
2618  s_fraction[2] = 0.0;
2619  // Set local node number
2620  n = i0 + n_p * i1;
2621  break;
2622 
2623  case 5:
2624  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2625  s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
2626  s_fraction[2] = 1.0;
2627  // Set local node number
2628  n = i0 + n_p * i1 + n_p * n_p * (n_p - 1);
2629  break;
2630  }
2631 
2632 
2633  // We have to check if the hi and lo directions along the
2634  // face are inverted or not
2635  Vector<double> s_in_neighb(3);
2636  for (unsigned i = 0; i < 3; i++)
2637  {
2638  s[i] = -1.0 + 2.0 * s_fraction[i];
2639  s_in_neighb[i] =
2640  s_lo_neigh[i] +
2641  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
2642  }
2643 
2644 
2645  // Loop over timesteps
2646  for (unsigned t = 0; t < n_time; t++)
2647  {
2648  // Get the nodal position from neighbour element
2649  Vector<double> x_in_neighb(3);
2650  neigh_pt->object_pt()->interpolated_x(
2651  t, s_in_neighb, x_in_neighb);
2652 
2653  // Check error only if the node is NOT periodic
2654  if (is_periodic == false)
2655  {
2656  // Check error
2657  for (unsigned i = 0; i < 3; i++)
2658  {
2659  // Find the spatial error
2660  double err =
2661  std::fabs(node_pt(n)->x(t, i) - x_in_neighb[i]);
2662 
2663  // If it's bigger than our tolerance, say so
2664  if (err > 1e-9)
2665  {
2666  oomph_info << "errx[" << i << "], t x, x_neigh: " << err
2667  << " " << t << " " << node_pt(n)->x(t, i)
2668  << " " << x_in_neighb[i] << std::endl;
2669  oomph_info << "at " << node_pt(n)->x(0) << " "
2670  << node_pt(n)->x(1) << " " << node_pt(n)->x(2)
2671  << " " << std::endl;
2672  }
2673 
2674  // If it's bigger than the previous max error, it is the
2675  // new max error!
2676  if (err > max_error_x[i])
2677  {
2678  max_error_x[i] = err;
2679  }
2680  }
2681  }
2682 
2683  // Get the values from neighbour element. Note: # of values
2684  // gets set by routine (because in general we don't know
2685  // how many interpolated values a certain element has
2686  Vector<double> values_in_neighb;
2687  neigh_pt->object_pt()->get_interpolated_values(
2688  t, s_in_neighb, values_in_neighb);
2689 
2690  // Get the values in current element.
2691  Vector<double> values;
2692  get_interpolated_values(t, s, values);
2693 
2694  // Now figure out how many continuously interpolated
2695  // values there are
2696  unsigned num_val =
2697  neigh_pt->object_pt()->ncont_interpolated_values();
2698 
2699  // Check error
2700  for (unsigned ival = 0; ival < num_val; ival++)
2701  {
2702  double err = std::fabs(values[ival] - values_in_neighb[ival]);
2703 
2704  if (err > 1.0e-10)
2705  {
2706  oomph_info << node_pt(n)->x(0) << " " << node_pt(n)->x(1)
2707  << " " << node_pt(n)->x(2) << " \n# "
2708  << "erru (S)" << err << " " << ival << " " << n
2709  << " " << values[ival] << " "
2710  << values_in_neighb[ival] << std::endl;
2711  // error_file<<"ZONE"<<std::endl
2712  // << node_pt(n)->x(0) << " "
2713  // << node_pt(n)->x(1) << " "
2714  // << node_pt(n)->x(2) << std::endl;
2715  }
2716 
2717  if (err > max_error_val)
2718  {
2719  max_error_val = err;
2720  }
2721  }
2722  }
2723  }
2724  }
2725  }
2726  }
2727  }
2728 
2729  max_error = max_error_x[0];
2730  if (max_error_x[1] > max_error) max_error = max_error_x[1];
2731  if (max_error_x[2] > max_error) max_error = max_error_x[2];
2732  if (max_error_val > max_error) max_error = max_error_val;
2733 
2734  if (max_error > 1e-9)
2735  {
2736  oomph_info << "\n#------------------------------------ \n#Max error ";
2737  oomph_info << max_error_x[0] << " " << max_error_x[1] << " "
2738  << max_error_x[2] << " " << max_error_val << std::endl;
2739  oomph_info << "#------------------------------------ \n " << std::endl;
2740  }
2741 
2742  // error_file.close();
2743  }
2744 
2745 
2746  //==================================================================
2747  /// Determine vector of solid (positional) boundary conditions along
2748  /// the element's boundary (or vertex) bound (S/W/N/E/SW/SE/NW/NE).
2749  ///
2750  /// This function assumes that the same boundary condition is applied
2751  /// along the entire length of an element's edge (of course, the
2752  /// vertices combine the boundary conditions of their two adjacent edges
2753  /// in the most restrictive combination. Hence, if we're at a vertex,
2754  /// we apply the most restrictive boundary condition of the
2755  /// two adjacent edges. If we're on an edge (in its proper interior),
2756  /// we apply the least restrictive boundary condition of all nodes
2757  /// along the edge.
2758  ///
2759  /// Usual convention:
2760  /// - solid_bound_cons[i]=0 if displacement in coordinate direction i
2761  /// on this boundary is free.
2762  /// - solid_bound_cons[i]=1 if it's pinned.
2763  //==================================================================
2764  void RefineableSolidQElement<3>::get_solid_bcs(
2765  int bound, Vector<int>& solid_bound_cons) const
2766  {
2767  using namespace OcTreeNames;
2768 
2769  // Spatial dimension of all nodes
2770  unsigned n_dim = this->nodal_dimension();
2771  // Set up temporary vectors to hold edge boundary conditions
2772  Vector<int> bound_cons1(n_dim), bound_cons2(n_dim);
2773  Vector<int> bound_cons3(n_dim);
2774 
2775  Vector<int> vect1(3), vect2(3), vect3(3);
2776  Vector<int> vect_elem;
2777  Vector<int> notzero;
2778  int n = 0;
2779 
2780  vect_elem = OcTree::Direction_to_vector[bound];
2781 
2782  // Just to see if bound is a face, an edge, or a vertex, n stores
2783  // the number of non-zero values in the vector reprensenting the bound
2784  // and the vector notzero stores the position of these values
2785  for (int i = 0; i < 3; i++)
2786  {
2787  if (vect_elem[i] != 0)
2788  {
2789  n++;
2790  notzero.push_back(i);
2791  }
2792  }
2793 
2794  switch (n)
2795  {
2796  // If there is only one non-zero value, bound is a face
2797  case 1:
2798  get_face_solid_bcs(bound, solid_bound_cons);
2799  break;
2800 
2801  // If there are two non-zero values, bound is an edge
2802  case 2:
2803 
2804  for (unsigned i = 0; i < 3; i++)
2805  {
2806  vect1[i] = 0;
2807  vect2[i] = 0;
2808  }
2809  // vect1 and vect2 are the vector of the two faces adjacent to bound
2810  vect1[notzero[0]] = vect_elem[notzero[0]];
2811  vect2[notzero[1]] = vect_elem[notzero[1]];
2812 
2813  get_face_solid_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
2814  get_face_solid_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
2815  // get the most restrictive bc
2816  for (unsigned k = 0; k < n_dim; k++)
2817  {
2818  solid_bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
2819  }
2820  break;
2821 
2822  // If there are three non-zero value, bound is a vertex
2823  case 3:
2824 
2825  for (unsigned i = 0; i < 3; i++)
2826  {
2827  vect1[i] = 0;
2828  vect2[i] = 0;
2829  vect3[i] = 0;
2830  }
2831  // vectors to the three adjacent faces of the vertex
2832  vect1[0] = vect_elem[0];
2833  vect2[1] = vect_elem[1];
2834  vect3[2] = vect_elem[2];
2835 
2836  get_face_solid_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
2837  get_face_solid_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
2838  get_face_solid_bcs(OcTree::Vector_to_direction[vect3], bound_cons3);
2839 
2840 
2841  // set the bcs to the most restrictive ones
2842  for (unsigned k = 0; k < n_dim; k++)
2843  {
2844  solid_bound_cons[k] =
2845  (bound_cons1[k] || bound_cons2[k] || bound_cons3[k]);
2846  }
2847  break;
2848 
2849  default:
2850  throw OomphLibError("Make sure you are not giving OMEGA as bound",
2851  OOMPH_CURRENT_FUNCTION,
2852  OOMPH_EXCEPTION_LOCATION);
2853  }
2854  }
2855 
2856  //==================================================================
2857  /// Determine Vector of solid (positional) boundary conditions along
2858  /// the element's edge (S/N/W/E) -- BC is the least restrictive combination
2859  /// of all the nodes on this edge
2860  ///
2861  /// Usual convention:
2862  /// - solid_bound_cons[i]=0 if displacement in coordinate direction i
2863  /// on this boundary is free
2864  /// - solid_bound_cons[i]=1 if it's pinned
2865  //==================================================================
2866  void RefineableSolidQElement<3>::get_face_solid_bcs(
2867  const int& face, Vector<int>& solid_bound_cons) const
2868  {
2869  using namespace OcTreeNames;
2870 
2871  // Number of nodes along 1D edge
2872  unsigned n_p = nnode_1d();
2873  // the four corner nodes on the boundary
2874  unsigned node1, node2, node3, node4;
2875 
2876  // Set the four corner nodes for the face
2877  switch (face)
2878  {
2879  case U:
2880  node1 = n_p * n_p * n_p - 1;
2881  node2 = n_p * n_p - 1;
2882  node3 = n_p * (n_p - 1);
2883  node4 = n_p * (n_p * n_p - 1);
2884  break;
2885 
2886  case D:
2887  node1 = 0;
2888  node2 = n_p - 1;
2889  node3 = (n_p * n_p + 1) * (n_p - 1);
2890  node4 = n_p * n_p * (n_p - 1);
2891  break;
2892 
2893  case R:
2894  node1 = n_p - 1;
2895  node2 = (n_p * n_p + 1) * (n_p - 1);
2896  node3 = n_p * n_p * n_p - 1;
2897  node4 = n_p * n_p - 1;
2898  break;
2899 
2900  case L:
2901  node1 = 0;
2902  node2 = n_p * (n_p - 1);
2903  node3 = n_p * (n_p * n_p - 1);
2904  node4 = n_p * n_p * (n_p - 1);
2905  break;
2906 
2907  case B:
2908  node1 = 0;
2909  node2 = n_p - 1;
2910  node3 = n_p * n_p - 1;
2911  node4 = n_p * (n_p - 1);
2912  break;
2913 
2914  case F:
2915  node1 = n_p * n_p * n_p - 1;
2916  node2 = n_p * (n_p * n_p - 1);
2917  node3 = n_p * n_p * (n_p - 1);
2918  node4 = (n_p - 1) * (n_p * n_p + 1);
2919  break;
2920 
2921  default:
2922  std::ostringstream error_stream;
2923  error_stream << "Wrong edge " << face << " passed\n";
2924 
2925  throw OomphLibError(
2926  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2927  }
2928 
2929  // Cast to solid nodes
2930  SolidNode* solid_node1_pt = dynamic_cast<SolidNode*>(node_pt(node1));
2931  SolidNode* solid_node2_pt = dynamic_cast<SolidNode*>(node_pt(node2));
2932  SolidNode* solid_node3_pt = dynamic_cast<SolidNode*>(node_pt(node3));
2933  SolidNode* solid_node4_pt = dynamic_cast<SolidNode*>(node_pt(node4));
2934 
2935  //#ifdef PARANOID
2936  if (solid_node1_pt == 0)
2937  {
2938  throw OomphLibError(
2939  "Corner node 1 cannot be cast to SolidNode --> something is wrong",
2940  OOMPH_CURRENT_FUNCTION,
2941  OOMPH_EXCEPTION_LOCATION);
2942  }
2943 
2944  if (solid_node2_pt == 0)
2945  {
2946  throw OomphLibError(
2947  "Corner node 2 cannot be cast to SolidNode --> something is wrong",
2948  OOMPH_CURRENT_FUNCTION,
2949  OOMPH_EXCEPTION_LOCATION);
2950  }
2951 
2952  if (solid_node3_pt == 0)
2953  {
2954  throw OomphLibError(
2955  "Corner node 3 cannot be cast to SolidNode --> something is wrong",
2956  OOMPH_CURRENT_FUNCTION,
2957  OOMPH_EXCEPTION_LOCATION);
2958  }
2959 
2960  if (solid_node4_pt == 0)
2961  {
2962  throw OomphLibError(
2963  "Corner node 4 cannot be cast to SolidNode --> something is wrong",
2964  OOMPH_CURRENT_FUNCTION,
2965  OOMPH_EXCEPTION_LOCATION);
2966  }
2967 
2968  //#endif
2969 
2970  // Number of coordinates
2971  unsigned n_dim = this->nodal_dimension();
2972 
2973  // Loop over all directions, the least restrictive value is
2974  // the multiplication of the boundary conditions at the 4 nodes
2975  // Assuming that free is always zero and pinned is one
2976  for (unsigned k = 0; k < n_dim; k++)
2977  {
2978  solid_bound_cons[k] = solid_node1_pt->position_is_pinned(k) *
2979  solid_node2_pt->position_is_pinned(k) *
2980  solid_node3_pt->position_is_pinned(k) *
2981  solid_node4_pt->position_is_pinned(k);
2982  }
2983  }
2984 
2985 
2986  //========================================================================
2987  /// Static matrix for coincidence between son nodal points and
2988  /// father boundaries
2989  ///
2990  //========================================================================
2991  std::map<unsigned, DenseMatrix<int>> RefineableQElement<3>::Father_bound;
2992 
2993 
2994 } // namespace oomph
algebraic_elements.h
e
e
Definition: cfortran.h:571
s
static char t char * s
Definition: cfortran.h:568
i
cstr elem_len * i
Definition: cfortran.h:603
t
char t
Definition: cfortran.h:568
oomph::AlgebraicElementBase
////////////////////////////////////////////////////////////////////
Definition: algebraic_elements.h:506
oomph::AlgebraicElementBase::setup_algebraic_node_update
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...
Definition: algebraic_elements.cc:48
oomph::Data::pin
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::Data::set_value
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:271
oomph::Data::nvalue
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
oomph::ElementWithMovingNodes
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: element_with_moving_nodes.h:65
oomph::ElementWithMovingNodes::enable_bypass_fill_in_jacobian_from_geometric_data
void enable_bypass_fill_in_jacobian_from_geometric_data()
Bypass the call to fill_in_jacobian_from_geometric_data.
Definition: element_with_moving_nodes.h:274
oomph::ElementWithMovingNodes::method_for_shape_derivs
int & method_for_shape_derivs()
Access to method (enumerated flag) for determination of shape derivs.
Definition: element_with_moving_nodes.h:268
oomph::ElementWithMovingNodes::enable_always_evaluate_dresidual_dnodal_coordinates_by_fd
void enable_always_evaluate_dresidual_dnodal_coordinates_by_fd()
Insist that shape derivatives are always evaluated by fd (using FiniteElement::get_dresidual_dnodal_c...
Definition: element_with_moving_nodes.h:194
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::get_node_at_local_coordinate
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3882
oomph::FiniteElement::s_min
virtual double s_min() const
Min value of local coordinate.
Definition: elements.h:2793
oomph::FiniteElement::interpolated_x
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:3962
oomph::FiniteElement::s_max
virtual double s_max() const
Max. value of local coordinate.
Definition: elements.h:2803
oomph::FiniteElement::get_x
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:1885
oomph::FiniteElement::macro_elem_pt
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Definition: elements.h:1878
oomph::HangInfo
Class that contains data for hanging nodes.
Definition: nodes.h:742
oomph::HangInfo::set_master_node_pt
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:1474
oomph::MacroElementNodeUpdateElementBase
////////////////////////////////////////////////////////////////////
Definition: macro_element_node_update_element.h:197
oomph::MacroElementNodeUpdateElementBase::geom_object_pt
Vector< GeomObject * > & geom_object_pt()
Vector of (pointers to) geometric objects involved in node update function.
Definition: macro_element_node_update_element.h:227
oomph::MacroElementNodeUpdateElementBase::set_node_update_info
virtual void set_node_update_info(const Vector< GeomObject * > &geom_object_pt)=0
Set node update information: Pass the vector of (pointers to) the geometric objects that affect the n...
oomph::Mesh
A general mesh class.
Definition: mesh.h:67
oomph::Mesh::add_boundary_node
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:243
oomph::Mesh::add_node_pt
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
oomph::Mesh::boundary_coordinate_exists
bool boundary_coordinate_exists(const unsigned &i) const
Indicate whether the i-th boundary has an intrinsic coordinate.
Definition: mesh.h:565
oomph::Node
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
oomph::Node::hanging_pt
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
oomph::Node::ndim
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
oomph::Node::set_coordinates_on_boundary
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:2394
oomph::Node::set_hanging_pt
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:2068
oomph::Node::make_periodic
virtual void make_periodic(Node *const &node_pt)
Make the node periodic by copying the values from node_pt. Note that the coordinates will always rema...
Definition: nodes.cc:2257
oomph::Node::nposition_type
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:1016
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
oomph::OcTree
OcTree class: Recursively defined, generalised octree.
Definition: octree.h:114
oomph::OcTree::faces_of_common_edge
static Vector< int > faces_of_common_edge(const int &edge)
Function that, given an edge, returns the two faces on which it.
Definition: octree.cc:268
oomph::OcTree::Direct_string
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: octree.h:329
oomph::OcTree::gteq_face_neighbour
OcTree * gteq_face_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_sw, Vector< double > &s_ne, int &face, int &diff_level, bool &in_neighbouring_tree) const
Find (pointer to) ‘greater-or-equal-sized face neighbour’ in given direction (L/R/U/D/F/B)....
Definition: octree.cc:3373
oomph::OcTree::gteq_true_edge_neighbour
OcTree * gteq_true_edge_neighbour(const int &direction, const unsigned &i_root_edge_neighbour, unsigned &nroot_edge_neighbour, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level) const
Find (pointer to) ‘greater-or-equal-sized true edge neighbour’ in the given direction (LB,...
Definition: octree.cc:3618
oomph::OcTree::Direction_to_vector
static Vector< Vector< int > > Direction_to_vector
For each direction, i.e. a son_type (vertex), a face or an edge, this defines a vector that indicates...
Definition: octree.h:353
oomph::OcTree::Vector_to_direction
static std::map< Vector< int >, int > Vector_to_direction
Each vector representing a direction can be translated into a direction, either a son type (vertex),...
Definition: octree.h:348
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition: oomph_definitions.h:222
oomph::OomphLibWarning
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Definition: oomph_definitions.h:267
oomph::QElementBase::s_macro_ll
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:186
oomph::QElementBase::s_macro_ur
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:202
oomph::RefineableElement
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
Definition: refineable_elements.h:97
oomph::RefineableElement::nodes_built
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
Definition: refineable_elements.h:389
oomph::RefineableElement::ncont_interpolated_values
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
oomph::RefineableElement::ninterpolating_node
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...
Definition: refineable_elements.h:466
oomph::RefineableElement::interpolating_basis
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...
Definition: refineable_elements.h:481
oomph::RefineableElement::interpolating_node_pt
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....
Definition: refineable_elements.h:437
oomph::RefineableElement::get_interpolating_node_at_local_coordinate
virtual Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &value_id)
Return a pointer to the node that interpolates the value-id-th unknown at local coordinate s....
Definition: refineable_elements.h:456
oomph::RefineableElement::get_interpolated_values
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...
Definition: refineable_elements.h:417
oomph::RefineableQElement< 3 >
Refineable version of QElement<3,NNODE_1D>.
Definition: refineable_brick_element.h:68
oomph::RefineableQElement< 3 >::interpolated_zeta_on_face
void interpolated_zeta_on_face(const unsigned &boundary, const int &face, const Vector< double > &s, Vector< double > &zeta)
Return the value of the intrinsic boundary coordinate interpolated along the face.
Definition: refineable_brick_element.cc:582
oomph::RefineableQElement< 3 >::get_bcs
void get_bcs(int bound, Vector< int > &bound_cons) const
Determine Vector of boundary conditions along the element's boundary (or vertex) bound (S/W/N/E/SW/SE...
Definition: refineable_brick_element.cc:238
oomph::RefineableQElement< 3 >::get_boundaries
void get_boundaries(const int &edge, std::set< unsigned > &boundaries) const
Given an element edge/vertex, return a Vector which contains all the (mesh-)boundary numbers that thi...
Definition: refineable_brick_element.cc:424
oomph::RefineableQElement
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
oomph::RefineableSolidQElement< 3 >
Refineable version of Solid brick elements.
Definition: refineable_brick_element.h:239
oomph::RefineableSolidQElement< 3 >::get_solid_bcs
void get_solid_bcs(int bound, Vector< int > &solid_bound_cons) const
Determine vector of solid (positional) boundary conditions along edge (or on vertex) bound (S/W/N/E/S...
Definition: refineable_brick_element.cc:2764
oomph::RefineableSolidQElement
A class that is used to template the solid refineable Q elements by dimension. It's really nothing mo...
Definition: Qelements.h:2286
oomph::Shape
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
oomph::SolidNode
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
oomph::SolidNode::pin_position
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
oomph::SolidNode::position_is_pinned
bool position_is_pinned(const unsigned &i)
Test whether the i-th coordinate is pinned, 0: false; 1: true.
Definition: nodes.h:1803
oomph::TimeStepper
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
oomph::TreeRoot::is_neighbour_periodic
bool is_neighbour_periodic(const int &direction)
Return whether the neighbour in the particular direction is periodic.
Definition: tree.h:364
oomph::Tree::object_pt
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:88
oomph::Tree::son_type
int son_type() const
Return son type.
Definition: tree.h:214
oomph::Tree::OMEGA
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:262
oomph::Tree::root_pt
TreeRoot *& root_pt()
Return pointer to root of the tree.
Definition: tree.h:141
macro_element_node_update_element.h
oomph::CumulativeTimings::start
void start(const unsigned &i)
(Re-)start i-th timer
Definition: oomph_utilities.cc:243
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
Definition: oomph_definitions.cc:286
oomph::Global_unsigned::Number
unsigned Number
The unsigned.
Definition: oomph_definitions.cc:269
oomph::OneDimLagrange::shape
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:564
oomph::TecplotNames::colour
Vector< std::string > colour
Tecplot colours.
Definition: oomph_utilities.cc:1159
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30
oomph::oomph_info
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...
Definition: oomph_definitions.cc:319
refineable_brick_element.h