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