brick_from_tet_mesh.template.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 #ifndef OOMPH_BRICK_FROM_TET_MESH_TEMPLATE_CC
31 #define OOMPH_BRICK_FROM_TET_MESH_TEMPLATE_CC
32 
33 
35 
36 namespace oomph
37 {
38 
39 
40  //=======================================================================
41  /// Build fct: Pass pointer to existing tet mesh and timestepper
42  /// Specialisation for XdaTetMesh<TElement<3,3> >
43  //=======================================================================
44  template<class ELEMENT>
46  XdaTetMesh<TElement<3,3> >* tet_mesh_pt,
47  TimeStepper* time_stepper_pt)
48  {
49  // Mesh can only be built with 3D Qelements.
50  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3, 3);
51 
52  // Figure out if the tet mesh is a solid mesh
53  bool tet_mesh_is_solid_mesh=false;
54  if (dynamic_cast<SolidFiniteElement*>(tet_mesh_pt->element_pt(0))!=0)
55  {
56  tet_mesh_is_solid_mesh=true;
57  }
58 
59  // Setup lookup scheme for local coordinates on triangular faces.
60  // The local coordinates identify the points on the triangular
61  // FaceElements on which we place the bottom layer of the
62  // brick nodes.
63  Vector<Vector<double> > s_face(19);
64  for (unsigned i=0;i<19;i++)
65  {
66  s_face[i].resize(2);
67 
68  switch (i)
69  {
70 
71  // Vertex nodes
72 
73  case 0:
74  s_face[i][0]=1.0;
75  s_face[i][1]=0.0;
76  break;
77 
78  case 1:
79  s_face[i][0]=0.0;
80  s_face[i][1]=1.0;
81  break;
82 
83  case 2:
84  s_face[i][0]=0.0;
85  s_face[i][1]=0.0;
86  break;
87 
88  // Midside nodes
89 
90  case 3:
91  s_face[i][0]=0.5;
92  s_face[i][1]=0.5;
93  break;
94 
95  case 4:
96  s_face[i][0]=0.0;
97  s_face[i][1]=0.5;
98  break;
99 
100  case 5:
101  s_face[i][0]=0.5;
102  s_face[i][1]=0.0;
103  break;
104 
105 
106  // Quarter side nodes
107 
108  case 6:
109  s_face[i][0]=0.75;
110  s_face[i][1]=0.25;
111  break;
112 
113  case 7:
114  s_face[i][0]=0.25;
115  s_face[i][1]=0.75;
116  break;
117 
118  case 8:
119  s_face[i][0]=0.0;
120  s_face[i][1]=0.75;
121  break;
122 
123  case 9:
124  s_face[i][0]=0.0;
125  s_face[i][1]=0.25;
126  break;
127 
128  case 10:
129  s_face[i][0]=0.25;
130  s_face[i][1]=0.0;
131  break;
132 
133  case 11:
134  s_face[i][0]=0.75;
135  s_face[i][1]=0.0;
136  break;
137 
138  // Central node
139 
140  case 12:
141  s_face[i][0]=1.0/3.0;
142  s_face[i][1]=1.0/3.0;
143  break;
144 
145 
146  // Vertical internal midside nodes connecting 2 and 3
147 
148  case 13:
149  s_face[i][0]=5.0/24.0;
150  s_face[i][1]=5.0/24.0;
151  break;
152 
153  case 14:
154  s_face[i][0]=5.0/12.0;
155  s_face[i][1]=5.0/12.0;
156  break;
157 
158  // Internal midside nodes connecting nodes 0 and 4
159 
160  case 15:
161  s_face[i][1]=5.0/24.0;
162  s_face[i][0]=7.0/12.0; // 1.0-2.0*5.0/24.0;
163  break;
164 
165  case 16:
166  s_face[i][1]=5.0/12.0;
167  s_face[i][0]=1.0/6.0; // 1.0-2.0*5.0/12.0;
168  break;
169 
170 
171  // Internal midside nodes connecting nodes 1 and 5
172 
173  case 17:
174  s_face[i][0]=5.0/24.0;
175  s_face[i][1]=7.0/12.0; // 1.0-2.0*5.0/24.0;
176  break;
177 
178  case 18:
179  s_face[i][0]=5.0/12.0;
180  s_face[i][1]=1.0/6.0; //1.0-2.0*5.0/12.0;
181  break;
182 
183  }
184  }
185 
186  // Set number of boundaries
187  unsigned nb=tet_mesh_pt->nboundary();
188  set_nboundary(nb);
189 
190  // Get ready for boundary lookup scheme
191  Boundary_element_pt.resize(nb);
192  Face_index_at_boundary.resize(nb);
193 
194  // Maps to check which nodes have already been done
195 
196  // Map that stores the new brick node corresponding to an existing tet node
197  std::map<Node*,Node*> tet_node_node_pt;
198 
199  // Map that stores node on an edge between two brick nodes
200  std::map<Edge,Node*> brick_edge_node_pt;
201 
202  // Map that stores node on face spanned by three tet nodes
203  std::map<TFace,Node*> tet_face_node_pt;
204 
205  // Create the four Dummy bricks:
206  //------------------------------
207  Vector<DummyBrickElement*> dummy_q_el_pt(4);
208  for (unsigned e=0;e<4;e++)
209  {
210  dummy_q_el_pt[e]=new DummyBrickElement;
211  for (unsigned j=0;j<8;j++)
212  {
213  dummy_q_el_pt[e]->construct_node(j);
214  }
215  }
216 
217  // Loop over the elements in the tet mesh
218  unsigned n_el_tet=tet_mesh_pt->nelement();
219  for (unsigned e_tet=0; e_tet<n_el_tet; e_tet++)
220  {
221  // Cast to ten-noded tet
222  TElement<3,3>* tet_el_pt=dynamic_cast<TElement<3,3>*>(
223  tet_mesh_pt->element_pt(e_tet));
224 
225 #ifdef PARANOID
226  if (tet_el_pt==0)
227  {
228  std::ostringstream error_stream;
229  error_stream
230  << "BrickFromTetMesh can only built from tet mesh containing\n"
231  << "ten-noded tets.\n";
232  throw OomphLibError(
233  error_stream.str(),
234  OOMPH_CURRENT_FUNCTION,
235  OOMPH_EXCEPTION_LOCATION);
236  }
237 #endif
238 
239  // Storage for the centroid node for this tet
240  Node* centroid_node_pt=0;
241 
242  // Internal mid brick-face nodes
243  Node* top_mid_face_node0_pt=0;
244  Node* right_mid_face_node0_pt=0;
245  Node* back_mid_face_node0_pt=0;
246 
247  Node* top_mid_face_node1_pt=0;
248  Node* right_mid_face_node1_pt=0;
249 
250  Node* top_mid_face_node2_pt=0;
251 
252  // Newly created brick elements
253  FiniteElement* brick_el0_pt=0;
254  FiniteElement* brick_el1_pt=0;
255  FiniteElement* brick_el2_pt=0;
256  FiniteElement* brick_el3_pt=0;
257 
258 
259  // First brick element is centred at node 0 of tet:
260  //-------------------------------------------------
261  {
262 
263  // Assign coordinates of dummy element
264  for (unsigned j=0;j<8;j++)
265  {
266  Node* nod_pt=dummy_q_el_pt[0]->node_pt(j);
267  Vector<double> s_tet(3);
268  Vector<double> x_tet(3);
269  switch (j)
270  {
271  case 0:
272  tet_el_pt->local_coordinate_of_node(0,s_tet);
273  nod_pt->set_value(0,s_tet[0]);
274  nod_pt->set_value(1,s_tet[1]);
275  nod_pt->set_value(2,s_tet[2]);
276  tet_el_pt->interpolated_x(s_tet,x_tet);
277  nod_pt->x(0)=x_tet[0];
278  nod_pt->x(1)=x_tet[1];
279  nod_pt->x(2)=x_tet[2];
280  break;
281  case 1:
282  tet_el_pt->local_coordinate_of_node(4,s_tet);
283  nod_pt->set_value(0,s_tet[0]);
284  nod_pt->set_value(1,s_tet[1]);
285  nod_pt->set_value(2,s_tet[2]);
286  tet_el_pt->interpolated_x(s_tet,x_tet);
287  nod_pt->x(0)=x_tet[0];
288  nod_pt->x(1)=x_tet[1];
289  nod_pt->x(2)=x_tet[2];
290  break;
291  case 2:
292  tet_el_pt->local_coordinate_of_node(6,s_tet);
293  nod_pt->set_value(0,s_tet[0]);
294  nod_pt->set_value(1,s_tet[1]);
295  nod_pt->set_value(2,s_tet[2]);
296  tet_el_pt->interpolated_x(s_tet,x_tet);
297  nod_pt->x(0)=x_tet[0];
298  nod_pt->x(1)=x_tet[1];
299  nod_pt->x(2)=x_tet[2];
300  break;
301  case 3:
302  // label 13 in initial sketch: Mid face node on face spanned by
303  // tet nodes 0,1,3
304  s_tet[0]=1.0/3.0;
305  s_tet[1]=1.0/3.0;
306  s_tet[2]=0.0;
307  nod_pt->set_value(0,s_tet[0]);
308  nod_pt->set_value(1,s_tet[1]);
309  nod_pt->set_value(2,s_tet[2]);
310  tet_el_pt->interpolated_x(s_tet,x_tet);
311  nod_pt->x(0)=x_tet[0];
312  nod_pt->x(1)=x_tet[1];
313  nod_pt->x(2)=x_tet[2];
314  break;
315  case 4:
316  tet_el_pt->local_coordinate_of_node(5,s_tet);
317  nod_pt->set_value(0,s_tet[0]);
318  nod_pt->set_value(1,s_tet[1]);
319  nod_pt->set_value(2,s_tet[2]);
320  tet_el_pt->interpolated_x(s_tet,x_tet);
321  nod_pt->x(0)=x_tet[0];
322  nod_pt->x(1)=x_tet[1];
323  nod_pt->x(2)=x_tet[2];
324  break;
325  case 5:
326  // label 11 in initial sketch: Mid face node on face spanned
327  // by tet nodes 0,1,2
328  s_tet[0]=1.0/3.0;
329  s_tet[1]=1.0/3.0;
330  s_tet[2]=1.0/3.0;
331  nod_pt->set_value(0,s_tet[0]);
332  nod_pt->set_value(1,s_tet[1]);
333  nod_pt->set_value(2,s_tet[2]);
334  tet_el_pt->interpolated_x(s_tet,x_tet);
335  nod_pt->x(0)=x_tet[0];
336  nod_pt->x(1)=x_tet[1];
337  nod_pt->x(2)=x_tet[2];
338  break;
339  case 6:
340  // label 12 in initial sketch: Mid face node on face
341  // spanned by tet nodes 0,2,3
342  s_tet[0]=1.0/3.0;
343  s_tet[1]=0.0;
344  s_tet[2]=1.0/3.0;
345  nod_pt->set_value(0,s_tet[0]);
346  nod_pt->set_value(1,s_tet[1]);
347  nod_pt->set_value(2,s_tet[2]);
348  tet_el_pt->interpolated_x(s_tet,x_tet);
349  nod_pt->x(0)=x_tet[0];
350  nod_pt->x(1)=x_tet[1];
351  nod_pt->x(2)=x_tet[2];
352  break;
353  case 7:
354  // label 14 in initial sketch: Centroid
355  s_tet[0]=0.25;
356  s_tet[1]=0.25;
357  s_tet[2]=0.25;
358  nod_pt->set_value(0,s_tet[0]);
359  nod_pt->set_value(1,s_tet[1]);
360  nod_pt->set_value(2,s_tet[2]);
361  tet_el_pt->interpolated_x(s_tet,x_tet);
362  nod_pt->x(0)=x_tet[0];
363  nod_pt->x(1)=x_tet[1];
364  nod_pt->x(2)=x_tet[2];
365  break;
366  }
367  }
368 
369 
370  // Create actual zeroth brick element
371  FiniteElement* el_pt=new ELEMENT;
372  brick_el0_pt=el_pt;
373  Element_pt.push_back(el_pt);
374 
375  TFace face0(tet_el_pt->node_pt(0),
376  tet_el_pt->node_pt(1),
377  tet_el_pt->node_pt(2));
378 
379  TFace face1(tet_el_pt->node_pt(0),
380  tet_el_pt->node_pt(2),
381  tet_el_pt->node_pt(3));
382 
383  TFace face2(tet_el_pt->node_pt(0),
384  tet_el_pt->node_pt(1),
385  tet_el_pt->node_pt(3));
386 
387 
388  // Tet vertex nodes along edges emanating from node 0 in brick
389  Vector<Vector<unsigned> > tet_edge_node(3);
390  tet_edge_node[0].resize(2);
391  tet_edge_node[0][0]=4;
392  tet_edge_node[0][1]=1;
393  tet_edge_node[1].resize(2);
394  tet_edge_node[1][0]=6;
395  tet_edge_node[1][1]=3;
396  tet_edge_node[2].resize(2);
397  tet_edge_node[2][0]=5;
398  tet_edge_node[2][1]=2;
399 
400  // Node number of tet vertex that node 0 in brick is centred on
401  unsigned central_tet_vertex=0;
402 
403  Node* tet_node_pt=0;
404  Node* old_node_pt=0;
405 
406  // Corner node
407  {
408  unsigned j=0;
409 
410  // Need new node?
411  tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
412  old_node_pt=tet_node_node_pt[tet_node_pt];
413  if (old_node_pt==0)
414  {
415  Node* new_node_pt=0;
416  if (tet_node_pt->is_on_boundary())
417  {
418  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
419  }
420  else
421  {
422  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
423  }
424  tet_node_node_pt[tet_node_pt]=new_node_pt;
425  Node_pt.push_back(new_node_pt);
426  Vector<double> s(3);
427  Vector<double> s_tet(3);
428  Vector<double> x_tet(3);
429  el_pt->local_coordinate_of_node(j,s);
430  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
431  tet_el_pt->interpolated_x(s_tet,x_tet);
432  new_node_pt->x(0)=x_tet[0];
433  new_node_pt->x(1)=x_tet[1];
434  new_node_pt->x(2)=x_tet[2];
435  }
436  // Node already exists
437  else
438  {
439  el_pt->node_pt(j)=old_node_pt;
440  }
441  }
442 
443 
444  // Brick vertex node coindides with mid-edge node on tet edge 0
445  {
446  unsigned j=2;
447 
448  // Need new node?
449  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
450  old_node_pt=tet_node_node_pt[tet_node_pt];
451  if (old_node_pt==0)
452  {
453  Node* new_node_pt=0;
454  if (tet_node_pt->is_on_boundary())
455  {
456  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
457  }
458  else
459  {
460  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
461  }
462  tet_node_node_pt[tet_node_pt]=new_node_pt;
463  Node_pt.push_back(new_node_pt);
464  Vector<double> s(3);
465  Vector<double> s_tet(3);
466  Vector<double> x_tet(3);
467  el_pt->local_coordinate_of_node(j,s);
468  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
469  tet_el_pt->interpolated_x(s_tet,x_tet);
470  new_node_pt->x(0)=x_tet[0];
471  new_node_pt->x(1)=x_tet[1];
472  new_node_pt->x(2)=x_tet[2];
473  }
474  // Node already exists
475  else
476  {
477  el_pt->node_pt(j)=old_node_pt;
478  }
479  }
480 
481 
482  // Brick vertex node coindides with mid vertex node of tet edge 1
483  {
484  unsigned j=6;
485 
486  // Need new node?
487  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
488  old_node_pt=tet_node_node_pt[tet_node_pt];
489  if (old_node_pt==0)
490  {
491  Node* new_node_pt=0;
492  if (tet_node_pt->is_on_boundary())
493  {
494  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
495  }
496  else
497  {
498  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
499  }
500  tet_node_node_pt[tet_node_pt]=new_node_pt;
501  Node_pt.push_back(new_node_pt);
502  Vector<double> s(3);
503  Vector<double> s_tet(3);
504  Vector<double> x_tet(3);
505  el_pt->local_coordinate_of_node(j,s);
506  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
507  tet_el_pt->interpolated_x(s_tet,x_tet);
508  new_node_pt->x(0)=x_tet[0];
509  new_node_pt->x(1)=x_tet[1];
510  new_node_pt->x(2)=x_tet[2];
511  }
512  // Node already exists
513  else
514  {
515  el_pt->node_pt(j)=old_node_pt;
516  }
517  }
518 
519 
520  // Brick vertex node coindides with mid-vertex node of tet edge 2
521  {
522  unsigned j=18;
523 
524  // Need new node?
525  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
526  old_node_pt=tet_node_node_pt[tet_node_pt];
527  if (old_node_pt==0)
528  {
529  Node* new_node_pt=0;
530  if (tet_node_pt->is_on_boundary())
531  {
532  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
533  }
534  else
535  {
536  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
537  }
538  tet_node_node_pt[tet_node_pt]=new_node_pt;
539  Node_pt.push_back(new_node_pt);
540  Vector<double> s(3);
541  Vector<double> s_tet(3);
542  Vector<double> x_tet(3);
543  el_pt->local_coordinate_of_node(j,s);
544  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
545  tet_el_pt->interpolated_x(s_tet,x_tet);
546  new_node_pt->x(0)=x_tet[0];
547  new_node_pt->x(1)=x_tet[1];
548  new_node_pt->x(2)=x_tet[2];
549  }
550  // Node already exists
551  else
552  {
553  el_pt->node_pt(j)=old_node_pt;
554  }
555  }
556 
557 
558 
559  // Brick vertex node in the middle of tet face0, spanned by
560  // tet vertices 0, 1, 2. Enumerated "11" in initial sketch.
561  {
562  unsigned j=20;
563 
564  // Need new node?
565  old_node_pt=tet_face_node_pt[face0];
566  if (old_node_pt==0)
567  {
568  Node* new_node_pt=0;
569  if (face0.is_boundary_face())
570  {
571  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
572  }
573  else
574  {
575  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
576  }
577  tet_face_node_pt[face0]=new_node_pt;
578  Node_pt.push_back(new_node_pt);
579  Vector<double> s(3);
580  Vector<double> s_tet(3);
581  Vector<double> x_tet(3);
582  el_pt->local_coordinate_of_node(j,s);
583  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
584  tet_el_pt->interpolated_x(s_tet,x_tet);
585  new_node_pt->x(0)=x_tet[0];
586  new_node_pt->x(1)=x_tet[1];
587  new_node_pt->x(2)=x_tet[2];
588  }
589  // Node already exists
590  else
591  {
592  el_pt->node_pt(j)=old_node_pt;
593  }
594  }
595 
596  // Brick vertex node in the middle of tet face1, spanned by
597  // tet vertices 0, 2, 3. Enumerated "12" in initial sketch.
598  {
599  unsigned j=24;
600 
601  // Need new node?
602  old_node_pt=tet_face_node_pt[face1];
603  if (old_node_pt==0)
604  {
605  Node* new_node_pt=0;
606  if (face1.is_boundary_face())
607  {
608  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
609  }
610  else
611  {
612  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
613  }
614  tet_face_node_pt[face1]=new_node_pt;
615  Node_pt.push_back(new_node_pt);
616  Vector<double> s(3);
617  Vector<double> s_tet(3);
618  Vector<double> x_tet(3);
619  el_pt->local_coordinate_of_node(j,s);
620  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
621  tet_el_pt->interpolated_x(s_tet,x_tet);
622  new_node_pt->x(0)=x_tet[0];
623  new_node_pt->x(1)=x_tet[1];
624  new_node_pt->x(2)=x_tet[2];
625  }
626  // Node already exists
627  else
628  {
629  el_pt->node_pt(j)=old_node_pt;
630  }
631  }
632 
633  // Brick vertex node in the middle of tet face2, spanned by
634  // tet vertices 0, 1, 3. Enumerated "13" in initial sketch.
635  {
636  unsigned j=8;
637 
638  // Need new node?
639  old_node_pt=tet_face_node_pt[face2];
640  if (old_node_pt==0)
641  {
642  Node* new_node_pt=0;
643  if (face2.is_boundary_face())
644  {
645  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
646  }
647  else
648  {
649  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
650  }
651  tet_face_node_pt[face2]=new_node_pt;
652  Node_pt.push_back(new_node_pt);
653  Vector<double> s(3);
654  Vector<double> s_tet(3);
655  Vector<double> x_tet(3);
656  el_pt->local_coordinate_of_node(j,s);
657  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
658  tet_el_pt->interpolated_x(s_tet,x_tet);
659  new_node_pt->x(0)=x_tet[0];
660  new_node_pt->x(1)=x_tet[1];
661  new_node_pt->x(2)=x_tet[2];
662  }
663  // Node already exists
664  else
665  {
666  el_pt->node_pt(j)=old_node_pt;
667  }
668  }
669 
670  // Brick vertex node in centroid of tet. Only built for first element.
671  // Enumerated "13" in initial sketch.
672  {
673  unsigned j=26;
674 
675  // Always new
676  {
677  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
678  centroid_node_pt=new_node_pt;
679  Node_pt.push_back(new_node_pt);
680  Vector<double> s(3);
681  Vector<double> s_tet(3);
682  Vector<double> x_tet(3);
683  el_pt->local_coordinate_of_node(j,s);
684  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
685  tet_el_pt->interpolated_x(s_tet,x_tet);
686  new_node_pt->x(0)=x_tet[0];
687  new_node_pt->x(1)=x_tet[1];
688  new_node_pt->x(2)=x_tet[2];
689  }
690  }
691 
692 
693  // Internal brick node -- always built
694  {
695  unsigned j=13;
696 
697  // Always new
698  {
699  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
700  Node_pt.push_back(new_node_pt);
701  Vector<double> s(3);
702  Vector<double> s_tet(3);
703  Vector<double> x_tet(3);
704  el_pt->local_coordinate_of_node(j,s);
705  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
706  tet_el_pt->interpolated_x(s_tet,x_tet);
707  new_node_pt->x(0)=x_tet[0];
708  new_node_pt->x(1)=x_tet[1];
709  new_node_pt->x(2)=x_tet[2];
710  }
711  }
712 
713  // Brick edge node between brick nodes 0 and 2
714  {
715  unsigned j=1;
716 
717  // Need new node?
718  Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
719  old_node_pt=brick_edge_node_pt[edge];
720  if (old_node_pt==0)
721  {
722  Node* new_node_pt=0;
723  if (edge.is_boundary_edge())
724  {
725  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
726  }
727  else
728  {
729  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
730  }
731  brick_edge_node_pt[edge]=new_node_pt;
732  Node_pt.push_back(new_node_pt);
733  Vector<double> s(3);
734  Vector<double> s_tet(3);
735  Vector<double> x_tet(3);
736  el_pt->local_coordinate_of_node(j,s);
737  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
738  tet_el_pt->interpolated_x(s_tet,x_tet);
739  new_node_pt->x(0)=x_tet[0];
740  new_node_pt->x(1)=x_tet[1];
741  new_node_pt->x(2)=x_tet[2];
742  }
743  // Node already exists
744  else
745  {
746  el_pt->node_pt(j)=old_node_pt;
747  }
748  }
749 
750 
751  // Brick edge node between brick nodes 0 and 6
752  {
753  unsigned j=3;
754 
755  // Need new node?
756  Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
757  old_node_pt=brick_edge_node_pt[edge];
758  if (old_node_pt==0)
759  {
760  Node* new_node_pt=0;
761  if (edge.is_boundary_edge())
762  {
763  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
764  }
765  else
766  {
767  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
768  }
769  brick_edge_node_pt[edge]=new_node_pt;
770  Node_pt.push_back(new_node_pt);
771  Vector<double> s(3);
772  Vector<double> s_tet(3);
773  Vector<double> x_tet(3);
774  el_pt->local_coordinate_of_node(j,s);
775  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
776  tet_el_pt->interpolated_x(s_tet,x_tet);
777  new_node_pt->x(0)=x_tet[0];
778  new_node_pt->x(1)=x_tet[1];
779  new_node_pt->x(2)=x_tet[2];
780  }
781  // Node already exists
782  else
783  {
784  el_pt->node_pt(j)=old_node_pt;
785  }
786  }
787 
788  // Brick edge node between brick nodes 2 and 8
789  {
790  unsigned j=5;
791 
792  // Need new node?
793  Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
794  old_node_pt=brick_edge_node_pt[edge];
795  if (old_node_pt==0)
796  {
797  Node* new_node_pt=0;
798  if (edge.is_boundary_edge())
799  {
800  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
801  }
802  else
803  {
804  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
805  }
806  brick_edge_node_pt[edge]=new_node_pt;
807  Node_pt.push_back(new_node_pt);
808  Vector<double> s(3);
809  Vector<double> s_tet(3);
810  Vector<double> x_tet(3);
811  el_pt->local_coordinate_of_node(j,s);
812  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
813  tet_el_pt->interpolated_x(s_tet,x_tet);
814  new_node_pt->x(0)=x_tet[0];
815  new_node_pt->x(1)=x_tet[1];
816  new_node_pt->x(2)=x_tet[2];
817  }
818  // Node already exists
819  else
820  {
821  el_pt->node_pt(j)=old_node_pt;
822  }
823  }
824 
825  // Brick edge node between brick nodes 6 and 8
826  {
827  unsigned j=7;
828 
829  // Need new node?
830  Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
831  old_node_pt=brick_edge_node_pt[edge];
832  if (old_node_pt==0)
833  {
834  Node* new_node_pt=0;
835  if (edge.is_boundary_edge())
836  {
837  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
838  }
839  else
840  {
841  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
842  }
843  brick_edge_node_pt[edge]=new_node_pt;
844  Node_pt.push_back(new_node_pt);
845  Vector<double> s(3);
846  Vector<double> s_tet(3);
847  Vector<double> x_tet(3);
848  el_pt->local_coordinate_of_node(j,s);
849  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
850  tet_el_pt->interpolated_x(s_tet,x_tet);
851  new_node_pt->x(0)=x_tet[0];
852  new_node_pt->x(1)=x_tet[1];
853  new_node_pt->x(2)=x_tet[2];
854  }
855  // Node already exists
856  else
857  {
858  el_pt->node_pt(j)=old_node_pt;
859  }
860  }
861 
862  // Brick edge node between brick nodes 18 and 20
863  {
864  unsigned j=19;
865 
866  // Need new node?
867  Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
868  old_node_pt=brick_edge_node_pt[edge];
869  if (old_node_pt==0)
870  {
871  Node* new_node_pt=0;
872  if (edge.is_boundary_edge())
873  {
874  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
875  }
876  else
877  {
878  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
879  }
880  brick_edge_node_pt[edge]=new_node_pt;
881  Node_pt.push_back(new_node_pt);
882  Vector<double> s(3);
883  Vector<double> s_tet(3);
884  Vector<double> x_tet(3);
885  el_pt->local_coordinate_of_node(j,s);
886  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
887  tet_el_pt->interpolated_x(s_tet,x_tet);
888  new_node_pt->x(0)=x_tet[0];
889  new_node_pt->x(1)=x_tet[1];
890  new_node_pt->x(2)=x_tet[2];
891  }
892  // Node already exists
893  else
894  {
895  el_pt->node_pt(j)=old_node_pt;
896  }
897  }
898 
899 
900  // Brick edge node between brick nodes 18 and 24
901  {
902  unsigned j=21;
903 
904  // Need new node?
905  Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
906  old_node_pt=brick_edge_node_pt[edge];
907  if (old_node_pt==0)
908  {
909  Node* new_node_pt=0;
910  if (edge.is_boundary_edge())
911  {
912  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
913  }
914  else
915  {
916  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
917  }
918  brick_edge_node_pt[edge]=new_node_pt;
919  Node_pt.push_back(new_node_pt);
920  Vector<double> s(3);
921  Vector<double> s_tet(3);
922  Vector<double> x_tet(3);
923  el_pt->local_coordinate_of_node(j,s);
924  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
925  tet_el_pt->interpolated_x(s_tet,x_tet);
926  new_node_pt->x(0)=x_tet[0];
927  new_node_pt->x(1)=x_tet[1];
928  new_node_pt->x(2)=x_tet[2];
929  }
930  // Node already exists
931  else
932  {
933  el_pt->node_pt(j)=old_node_pt;
934  }
935  }
936 
937  // Brick edge node between brick nodes 20 and 26
938  {
939  unsigned j=23;
940 
941  // Need new node?
942  Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
943  old_node_pt=brick_edge_node_pt[edge];
944  if (old_node_pt==0)
945  {
946  Node* new_node_pt=0;
947  if (edge.is_boundary_edge())
948  {
949  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
950  }
951  else
952  {
953  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
954  }
955  brick_edge_node_pt[edge]=new_node_pt;
956  Node_pt.push_back(new_node_pt);
957  Vector<double> s(3);
958  Vector<double> s_tet(3);
959  Vector<double> x_tet(3);
960  el_pt->local_coordinate_of_node(j,s);
961  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
962  tet_el_pt->interpolated_x(s_tet,x_tet);
963  new_node_pt->x(0)=x_tet[0];
964  new_node_pt->x(1)=x_tet[1];
965  new_node_pt->x(2)=x_tet[2];
966  }
967  // Node already exists
968  else
969  {
970  el_pt->node_pt(j)=old_node_pt;
971  }
972  }
973 
974 
975  // Brick edge node between brick nodes 24 and 26
976  {
977  unsigned j=25;
978 
979  // Need new node?
980  Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
981  old_node_pt=brick_edge_node_pt[edge];
982  if (old_node_pt==0)
983  {
984  Node* new_node_pt=0;
985  if (edge.is_boundary_edge())
986  {
987  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
988  }
989  else
990  {
991  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
992  }
993  brick_edge_node_pt[edge]=new_node_pt;
994  Node_pt.push_back(new_node_pt);
995  Vector<double> s(3);
996  Vector<double> s_tet(3);
997  Vector<double> x_tet(3);
998  el_pt->local_coordinate_of_node(j,s);
999  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1000  tet_el_pt->interpolated_x(s_tet,x_tet);
1001  new_node_pt->x(0)=x_tet[0];
1002  new_node_pt->x(1)=x_tet[1];
1003  new_node_pt->x(2)=x_tet[2];
1004  }
1005  // Node already exists
1006  else
1007  {
1008  el_pt->node_pt(j)=old_node_pt;
1009  }
1010  }
1011 
1012  // Brick edge node between brick nodes 0 and 18
1013  {
1014  unsigned j=9;
1015 
1016  // Need new node?
1017  Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
1018  old_node_pt=brick_edge_node_pt[edge];
1019  if (old_node_pt==0)
1020  {
1021  Node* new_node_pt=0;
1022  if (edge.is_boundary_edge())
1023  {
1024  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1025  }
1026  else
1027  {
1028  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1029  }
1030  brick_edge_node_pt[edge]=new_node_pt;
1031  Node_pt.push_back(new_node_pt);
1032  Vector<double> s(3);
1033  Vector<double> s_tet(3);
1034  Vector<double> x_tet(3);
1035  el_pt->local_coordinate_of_node(j,s);
1036  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1037  tet_el_pt->interpolated_x(s_tet,x_tet);
1038  new_node_pt->x(0)=x_tet[0];
1039  new_node_pt->x(1)=x_tet[1];
1040  new_node_pt->x(2)=x_tet[2];
1041  }
1042  // Node already exists
1043  else
1044  {
1045  el_pt->node_pt(j)=old_node_pt;
1046  }
1047  }
1048 
1049 
1050  // Brick edge node between brick nodes 2 and 20
1051  {
1052  unsigned j=11;
1053 
1054  // Need new node?
1055  Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
1056  old_node_pt=brick_edge_node_pt[edge];
1057  if (old_node_pt==0)
1058  {
1059  Node* new_node_pt=0;
1060  if (edge.is_boundary_edge())
1061  {
1062  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1063  }
1064  else
1065  {
1066  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1067  }
1068  brick_edge_node_pt[edge]=new_node_pt;
1069  Node_pt.push_back(new_node_pt);
1070  Vector<double> s(3);
1071  Vector<double> s_tet(3);
1072  Vector<double> x_tet(3);
1073  el_pt->local_coordinate_of_node(j,s);
1074  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1075  tet_el_pt->interpolated_x(s_tet,x_tet);
1076  new_node_pt->x(0)=x_tet[0];
1077  new_node_pt->x(1)=x_tet[1];
1078  new_node_pt->x(2)=x_tet[2];
1079  }
1080  // Node already exists
1081  else
1082  {
1083  el_pt->node_pt(j)=old_node_pt;
1084  }
1085  }
1086 
1087 
1088  // Brick edge node between brick nodes 6 and 24
1089  {
1090  unsigned j=15;
1091 
1092  // Need new node?
1093  Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
1094  old_node_pt=brick_edge_node_pt[edge];
1095  if (old_node_pt==0)
1096  {
1097  Node* new_node_pt=0;
1098  if (edge.is_boundary_edge())
1099  {
1100  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1101  }
1102  else
1103  {
1104  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1105  }
1106  brick_edge_node_pt[edge]=new_node_pt;
1107  Node_pt.push_back(new_node_pt);
1108  Vector<double> s(3);
1109  Vector<double> s_tet(3);
1110  Vector<double> x_tet(3);
1111  el_pt->local_coordinate_of_node(j,s);
1112  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1113  tet_el_pt->interpolated_x(s_tet,x_tet);
1114  new_node_pt->x(0)=x_tet[0];
1115  new_node_pt->x(1)=x_tet[1];
1116  new_node_pt->x(2)=x_tet[2];
1117  }
1118  // Node already exists
1119  else
1120  {
1121  el_pt->node_pt(j)=old_node_pt;
1122  }
1123  }
1124 
1125 
1126  // Brick edge node between brick nodes 8 and 26
1127  {
1128  unsigned j=17;
1129 
1130  // Need new node?
1131  Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
1132  old_node_pt=brick_edge_node_pt[edge];
1133  if (old_node_pt==0)
1134  {
1135  Node* new_node_pt=0;
1136  if (edge.is_boundary_edge())
1137  {
1138  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1139  }
1140  else
1141  {
1142  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1143  }
1144  brick_edge_node_pt[edge]=new_node_pt;
1145  Node_pt.push_back(new_node_pt);
1146  Vector<double> s(3);
1147  Vector<double> s_tet(3);
1148  Vector<double> x_tet(3);
1149  el_pt->local_coordinate_of_node(j,s);
1150  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1151  tet_el_pt->interpolated_x(s_tet,x_tet);
1152  new_node_pt->x(0)=x_tet[0];
1153  new_node_pt->x(1)=x_tet[1];
1154  new_node_pt->x(2)=x_tet[2];
1155  }
1156  // Node already exists
1157  else
1158  {
1159  el_pt->node_pt(j)=old_node_pt;
1160  }
1161  }
1162 
1163 
1164  // Mid brick-face node associated with face
1165  // spanned by mid-vertex nodes associated with tet edges 0 and 2
1166  {
1167  unsigned j=10;
1168 
1169  // Need new node?
1170  TFace face(tet_el_pt->node_pt( central_tet_vertex),
1171  tet_el_pt->node_pt(tet_edge_node[0][0]),
1172  tet_el_pt->node_pt(tet_edge_node[2][0]));
1173 
1174  old_node_pt=tet_face_node_pt[face];
1175  if (old_node_pt==0)
1176  {
1177  Node* new_node_pt=0;
1178  if (face.is_boundary_face())
1179  {
1180  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1181  }
1182  else
1183  {
1184  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1185  }
1186  tet_face_node_pt[face]=new_node_pt;
1187  Node_pt.push_back(new_node_pt);
1188  Vector<double> s(3);
1189  Vector<double> s_tet(3);
1190  Vector<double> x_tet(3);
1191  el_pt->local_coordinate_of_node(j,s);
1192  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1193  tet_el_pt->interpolated_x(s_tet,x_tet);
1194  new_node_pt->x(0)=x_tet[0];
1195  new_node_pt->x(1)=x_tet[1];
1196  new_node_pt->x(2)=x_tet[2];
1197  }
1198  // Node already exists
1199  else
1200  {
1201  el_pt->node_pt(j)=old_node_pt;
1202  }
1203  }
1204 
1205 
1206 
1207  // Mid brick-face node associated with face
1208  // spanned by mid-vertex nodes associated with tet edges 1 and 2
1209  {
1210  unsigned j=12;
1211 
1212  // Need new node?
1213  TFace face(tet_el_pt->node_pt(central_tet_vertex),
1214  tet_el_pt->node_pt(tet_edge_node[1][0]),
1215  tet_el_pt->node_pt(tet_edge_node[2][0]));
1216 
1217  old_node_pt=tet_face_node_pt[face];
1218  if (old_node_pt==0)
1219  {
1220  Node* new_node_pt=0;
1221  if (face.is_boundary_face())
1222  {
1223  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1224  }
1225  else
1226  {
1227  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1228  }
1229  tet_face_node_pt[face]=new_node_pt;
1230  Node_pt.push_back(new_node_pt);
1231  Vector<double> s(3);
1232  Vector<double> s_tet(3);
1233  Vector<double> x_tet(3);
1234  el_pt->local_coordinate_of_node(j,s);
1235  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1236  tet_el_pt->interpolated_x(s_tet,x_tet);
1237  new_node_pt->x(0)=x_tet[0];
1238  new_node_pt->x(1)=x_tet[1];
1239  new_node_pt->x(2)=x_tet[2];
1240  }
1241  // Node already exists
1242  else
1243  {
1244  el_pt->node_pt(j)=old_node_pt;
1245  }
1246  }
1247 
1248 
1249 
1250  // Mid brick-face node associated with face
1251  // spanned by mid-vertex nodes associated with tet edges 0 and 1
1252  {
1253  unsigned j=4;
1254 
1255  // Need new node?
1256  TFace face(tet_el_pt->node_pt(central_tet_vertex),
1257  tet_el_pt->node_pt(tet_edge_node[0][0]),
1258  tet_el_pt->node_pt(tet_edge_node[1][0]));
1259 
1260  old_node_pt=tet_face_node_pt[face];
1261  if (old_node_pt==0)
1262  {
1263  Node* new_node_pt=0;
1264  if (face.is_boundary_face())
1265  {
1266  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1267  }
1268  else
1269  {
1270  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1271  }
1272  tet_face_node_pt[face]=new_node_pt;
1273  Node_pt.push_back(new_node_pt);
1274  Vector<double> s(3);
1275  Vector<double> s_tet(3);
1276  Vector<double> x_tet(3);
1277  el_pt->local_coordinate_of_node(j,s);
1278  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1279  tet_el_pt->interpolated_x(s_tet,x_tet);
1280  new_node_pt->x(0)=x_tet[0];
1281  new_node_pt->x(1)=x_tet[1];
1282  new_node_pt->x(2)=x_tet[2];
1283  }
1284  // Node already exists
1285  else
1286  {
1287  el_pt->node_pt(j)=old_node_pt;
1288  }
1289  }
1290 
1291 
1292  // Top mid brick-face node -- only built by first element
1293  {
1294  unsigned j=22;
1295  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1296  Node_pt.push_back(new_node_pt);
1297  Vector<double> s(3);
1298  Vector<double> s_tet(3);
1299  Vector<double> x_tet(3);
1300  el_pt->local_coordinate_of_node(j,s);
1301  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1302  top_mid_face_node0_pt=new_node_pt;
1303  tet_el_pt->interpolated_x(s_tet,x_tet);
1304  new_node_pt->x(0)=x_tet[0];
1305  new_node_pt->x(1)=x_tet[1];
1306  new_node_pt->x(2)=x_tet[2];
1307  }
1308 
1309 
1310 
1311  // Right mid brick-face node -- only built by first element
1312  {
1313  unsigned j=14;
1314  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1315  Node_pt.push_back(new_node_pt);
1316  Vector<double> s(3);
1317  Vector<double> s_tet(3);
1318  Vector<double> x_tet(3);
1319  el_pt->local_coordinate_of_node(j,s);
1320  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1321  right_mid_face_node0_pt=new_node_pt;
1322  tet_el_pt->interpolated_x(s_tet,x_tet);
1323  new_node_pt->x(0)=x_tet[0];
1324  new_node_pt->x(1)=x_tet[1];
1325  new_node_pt->x(2)=x_tet[2];
1326  }
1327 
1328 
1329  // Back mid brick-face node -- only built by first element
1330  {
1331  unsigned j=16;
1332  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1333  Node_pt.push_back(new_node_pt);
1334  Vector<double> s(3);
1335  Vector<double> s_tet(3);
1336  Vector<double> x_tet(3);
1337  el_pt->local_coordinate_of_node(j,s);
1338  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1339  back_mid_face_node0_pt=new_node_pt;
1340  tet_el_pt->interpolated_x(s_tet,x_tet);
1341  new_node_pt->x(0)=x_tet[0];
1342  new_node_pt->x(1)=x_tet[1];
1343  new_node_pt->x(2)=x_tet[2];
1344  }
1345 
1346  }
1347 
1348 
1349  // Second brick element is centred at node 1 of tet:
1350  //--------------------------------------------------
1351  {
1352  // Assign coordinates of dummy element
1353  for (unsigned j=0;j<8;j++)
1354  {
1355  Node* nod_pt=dummy_q_el_pt[1]->node_pt(j);
1356  Vector<double> s_tet(3);
1357  Vector<double> x_tet(3);
1358  switch (j)
1359  {
1360  case 0:
1361  tet_el_pt->local_coordinate_of_node(1,s_tet);
1362  nod_pt->set_value(0,s_tet[0]);
1363  nod_pt->set_value(1,s_tet[1]);
1364  nod_pt->set_value(2,s_tet[2]);
1365  tet_el_pt->interpolated_x(s_tet,x_tet);
1366  nod_pt->x(0)=x_tet[0];
1367  nod_pt->x(1)=x_tet[1];
1368  nod_pt->x(2)=x_tet[2];
1369  break;
1370  case 1:
1371  tet_el_pt->local_coordinate_of_node(9,s_tet);
1372  nod_pt->set_value(0,s_tet[0]);
1373  nod_pt->set_value(1,s_tet[1]);
1374  nod_pt->set_value(2,s_tet[2]);
1375  tet_el_pt->interpolated_x(s_tet,x_tet);
1376  nod_pt->x(0)=x_tet[0];
1377  nod_pt->x(1)=x_tet[1];
1378  nod_pt->x(2)=x_tet[2];
1379  break;
1380  case 2:
1381  tet_el_pt->local_coordinate_of_node(4,s_tet);
1382  nod_pt->set_value(0,s_tet[0]);
1383  nod_pt->set_value(1,s_tet[1]);
1384  nod_pt->set_value(2,s_tet[2]);
1385  tet_el_pt->interpolated_x(s_tet,x_tet);
1386  nod_pt->x(0)=x_tet[0];
1387  nod_pt->x(1)=x_tet[1];
1388  nod_pt->x(2)=x_tet[2];
1389  break;
1390  case 3:
1391  // label 13 in initial sketch: Mid face node on face
1392  // spanned by tet nodes 0,1,3
1393  s_tet[0]=1.0/3.0;
1394  s_tet[1]=1.0/3.0;
1395  s_tet[2]=0.0;
1396  nod_pt->set_value(0,s_tet[0]);
1397  nod_pt->set_value(1,s_tet[1]);
1398  nod_pt->set_value(2,s_tet[2]);
1399  tet_el_pt->interpolated_x(s_tet,x_tet);
1400  nod_pt->x(0)=x_tet[0];
1401  nod_pt->x(1)=x_tet[1];
1402  nod_pt->x(2)=x_tet[2];
1403  break;
1404  case 4:
1405  tet_el_pt->local_coordinate_of_node(7,s_tet);
1406  nod_pt->set_value(0,s_tet[0]);
1407  nod_pt->set_value(1,s_tet[1]);
1408  nod_pt->set_value(2,s_tet[2]);
1409  tet_el_pt->interpolated_x(s_tet,x_tet);
1410  nod_pt->x(0)=x_tet[0];
1411  nod_pt->x(1)=x_tet[1];
1412  nod_pt->x(2)=x_tet[2];
1413  break;
1414  case 5:
1415  // label 10 in initial sketch: Mid face node on face
1416  // spanned by tet nodes 1,2,3
1417  s_tet[0]=0.0;
1418  s_tet[1]=1.0/3.0;
1419  s_tet[2]=1.0/3.0;
1420  nod_pt->set_value(0,s_tet[0]);
1421  nod_pt->set_value(1,s_tet[1]);
1422  nod_pt->set_value(2,s_tet[2]);
1423  tet_el_pt->interpolated_x(s_tet,x_tet);
1424  nod_pt->x(0)=x_tet[0];
1425  nod_pt->x(1)=x_tet[1];
1426  nod_pt->x(2)=x_tet[2];
1427  break;
1428  case 6:
1429  // label 11 in initial sketch: Mid face node on face
1430  // spanned by tet nodes 0,1,2
1431  s_tet[0]=1.0/3.0;
1432  s_tet[1]=1.0/3.0;
1433  s_tet[2]=1.0/3.0;
1434  nod_pt->set_value(0,s_tet[0]);
1435  nod_pt->set_value(1,s_tet[1]);
1436  nod_pt->set_value(2,s_tet[2]);
1437  tet_el_pt->interpolated_x(s_tet,x_tet);
1438  nod_pt->x(0)=x_tet[0];
1439  nod_pt->x(1)=x_tet[1];
1440  nod_pt->x(2)=x_tet[2];
1441  break;
1442  case 7:
1443  // label 14 in initial sketch: Centroid
1444  s_tet[0]=0.25;
1445  s_tet[1]=0.25;
1446  s_tet[2]=0.25;
1447  nod_pt->set_value(0,s_tet[0]);
1448  nod_pt->set_value(1,s_tet[1]);
1449  nod_pt->set_value(2,s_tet[2]);
1450  tet_el_pt->interpolated_x(s_tet,x_tet);
1451  nod_pt->x(0)=x_tet[0];
1452  nod_pt->x(1)=x_tet[1];
1453  nod_pt->x(2)=x_tet[2];
1454  break;
1455  }
1456  }
1457 
1458 
1459  // Create actual first brick element
1460  FiniteElement* el_pt=new ELEMENT;
1461  brick_el1_pt=el_pt;
1462  Element_pt.push_back(el_pt);
1463 
1464  TFace face0(tet_el_pt->node_pt(1),
1465  tet_el_pt->node_pt(3),
1466  tet_el_pt->node_pt(2));
1467 
1468  TFace face1(tet_el_pt->node_pt(1),
1469  tet_el_pt->node_pt(0),
1470  tet_el_pt->node_pt(2));
1471 
1472  TFace face2(tet_el_pt->node_pt(1),
1473  tet_el_pt->node_pt(0),
1474  tet_el_pt->node_pt(3));
1475 
1476  // Tet vertex nodes along edges emanating from node 0 in brick
1477  Vector<Vector<unsigned> > tet_edge_node(3);
1478  tet_edge_node[0].resize(2);
1479  tet_edge_node[0][0]=9;
1480  tet_edge_node[0][1]=3;
1481  tet_edge_node[1].resize(2);
1482  tet_edge_node[1][0]=4;
1483  tet_edge_node[1][1]=0;
1484  tet_edge_node[2].resize(2);
1485  tet_edge_node[2][0]=7;
1486  tet_edge_node[2][1]=2;
1487 
1488  // Node number of tet vertex that node 0 in brick is centred on
1489  unsigned central_tet_vertex=1;
1490 
1491  Node* tet_node_pt=0;
1492  Node* old_node_pt=0;
1493 
1494  // Corner node
1495  {
1496  unsigned j=0;
1497 
1498  // Need new node?
1499  tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
1500  old_node_pt=tet_node_node_pt[tet_node_pt];
1501  if (old_node_pt==0)
1502  {
1503  Node* new_node_pt=0;
1504  if (tet_node_pt->is_on_boundary())
1505  {
1506  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1507  }
1508  else
1509  {
1510  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1511  }
1512  tet_node_node_pt[tet_node_pt]=new_node_pt;
1513  Node_pt.push_back(new_node_pt);
1514  Vector<double> s(3);
1515  Vector<double> s_tet(3);
1516  Vector<double> x_tet(3);
1517  el_pt->local_coordinate_of_node(j,s);
1518  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1519  tet_el_pt->interpolated_x(s_tet,x_tet);
1520  new_node_pt->x(0)=x_tet[0];
1521  new_node_pt->x(1)=x_tet[1];
1522  new_node_pt->x(2)=x_tet[2];
1523  }
1524  // Node already exists
1525  else
1526  {
1527  el_pt->node_pt(j)=old_node_pt;
1528  }
1529  }
1530 
1531 
1532  // Brick vertex node coindides with mid-edge node on tet edge 0
1533  {
1534  unsigned j=2;
1535 
1536  // Need new node?
1537  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
1538  old_node_pt=tet_node_node_pt[tet_node_pt];
1539  if (old_node_pt==0)
1540  {
1541  Node* new_node_pt=0;
1542  if (tet_node_pt->is_on_boundary())
1543  {
1544  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1545  }
1546  else
1547  {
1548  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1549  }
1550  tet_node_node_pt[tet_node_pt]=new_node_pt;
1551  Node_pt.push_back(new_node_pt);
1552  Vector<double> s(3);
1553  Vector<double> s_tet(3);
1554  Vector<double> x_tet(3);
1555  el_pt->local_coordinate_of_node(j,s);
1556  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1557  tet_el_pt->interpolated_x(s_tet,x_tet);
1558  new_node_pt->x(0)=x_tet[0];
1559  new_node_pt->x(1)=x_tet[1];
1560  new_node_pt->x(2)=x_tet[2];
1561  }
1562  // Node already exists
1563  else
1564  {
1565  el_pt->node_pt(j)=old_node_pt;
1566  }
1567  }
1568 
1569 
1570  // Brick vertex node coindides with mid vertex node of tet edge 1
1571  {
1572  unsigned j=6;
1573 
1574  // Need new node?
1575  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
1576  old_node_pt=tet_node_node_pt[tet_node_pt];
1577  if (old_node_pt==0)
1578  {
1579  Node* new_node_pt=0;
1580  if (tet_node_pt->is_on_boundary())
1581  {
1582  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1583  }
1584  else
1585  {
1586  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1587  }
1588  tet_node_node_pt[tet_node_pt]=new_node_pt;
1589  Node_pt.push_back(new_node_pt);
1590  Vector<double> s(3);
1591  Vector<double> s_tet(3);
1592  Vector<double> x_tet(3);
1593  el_pt->local_coordinate_of_node(j,s);
1594  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1595  tet_el_pt->interpolated_x(s_tet,x_tet);
1596  new_node_pt->x(0)=x_tet[0];
1597  new_node_pt->x(1)=x_tet[1];
1598  new_node_pt->x(2)=x_tet[2];
1599  }
1600  // Node already exists
1601  else
1602  {
1603  el_pt->node_pt(j)=old_node_pt;
1604  }
1605  }
1606 
1607 
1608  // Brick vertex node coindides with mid-vertex node of tet edge 2
1609  {
1610  unsigned j=18;
1611 
1612  // Need new node?
1613  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
1614  old_node_pt=tet_node_node_pt[tet_node_pt];
1615  if (old_node_pt==0)
1616  {
1617  Node* new_node_pt=0;
1618  if (tet_node_pt->is_on_boundary())
1619  {
1620  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1621  }
1622  else
1623  {
1624  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1625  }
1626  tet_node_node_pt[tet_node_pt]=new_node_pt;
1627  Node_pt.push_back(new_node_pt);
1628  Vector<double> s(3);
1629  Vector<double> s_tet(3);
1630  Vector<double> x_tet(3);
1631  el_pt->local_coordinate_of_node(j,s);
1632  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1633  tet_el_pt->interpolated_x(s_tet,x_tet);
1634  new_node_pt->x(0)=x_tet[0];
1635  new_node_pt->x(1)=x_tet[1];
1636  new_node_pt->x(2)=x_tet[2];
1637  }
1638  // Node already exists
1639  else
1640  {
1641  el_pt->node_pt(j)=old_node_pt;
1642  }
1643  }
1644 
1645 
1646 
1647  // Brick vertex node in the middle of tet face0
1648  {
1649  unsigned j=20;
1650 
1651  // Need new node?
1652  old_node_pt=tet_face_node_pt[face0];
1653  if (old_node_pt==0)
1654  {
1655  Node* new_node_pt=0;
1656  if (face0.is_boundary_face())
1657  {
1658  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1659  }
1660  else
1661  {
1662  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1663  }
1664  tet_face_node_pt[face0]=new_node_pt;
1665  Node_pt.push_back(new_node_pt);
1666  Vector<double> s(3);
1667  Vector<double> s_tet(3);
1668  Vector<double> x_tet(3);
1669  el_pt->local_coordinate_of_node(j,s);
1670  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1671  tet_el_pt->interpolated_x(s_tet,x_tet);
1672  new_node_pt->x(0)=x_tet[0];
1673  new_node_pt->x(1)=x_tet[1];
1674  new_node_pt->x(2)=x_tet[2];
1675  }
1676  // Node already exists
1677  else
1678  {
1679  el_pt->node_pt(j)=old_node_pt;
1680  }
1681  }
1682 
1683  // Brick vertex node in the middle of tet face1
1684  {
1685  unsigned j=24;
1686 
1687  // Need new node?
1688  old_node_pt=tet_face_node_pt[face1];
1689  if (old_node_pt==0)
1690  {
1691  Node* new_node_pt=0;
1692  if (face1.is_boundary_face())
1693  {
1694  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1695  }
1696  else
1697  {
1698  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1699  }
1700  tet_face_node_pt[face1]=new_node_pt;
1701  Node_pt.push_back(new_node_pt);
1702  Vector<double> s(3);
1703  Vector<double> s_tet(3);
1704  Vector<double> x_tet(3);
1705  el_pt->local_coordinate_of_node(j,s);
1706  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1707  tet_el_pt->interpolated_x(s_tet,x_tet);
1708  new_node_pt->x(0)=x_tet[0];
1709  new_node_pt->x(1)=x_tet[1];
1710  new_node_pt->x(2)=x_tet[2];
1711  }
1712  // Node already exists
1713  else
1714  {
1715  el_pt->node_pt(j)=old_node_pt;
1716  }
1717  }
1718 
1719  // Brick vertex node in the middle of face2
1720  {
1721  unsigned j=8;
1722 
1723  // Need new node?
1724  old_node_pt=tet_face_node_pt[face2];
1725  if (old_node_pt==0)
1726  {
1727  Node* new_node_pt=0;
1728  if (face2.is_boundary_face())
1729  {
1730  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1731  }
1732  else
1733  {
1734  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1735  }
1736  tet_face_node_pt[face2]=new_node_pt;
1737  Node_pt.push_back(new_node_pt);
1738  Vector<double> s(3);
1739  Vector<double> s_tet(3);
1740  Vector<double> x_tet(3);
1741  el_pt->local_coordinate_of_node(j,s);
1742  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1743  tet_el_pt->interpolated_x(s_tet,x_tet);
1744  new_node_pt->x(0)=x_tet[0];
1745  new_node_pt->x(1)=x_tet[1];
1746  new_node_pt->x(2)=x_tet[2];
1747  }
1748  // Node already exists
1749  else
1750  {
1751  el_pt->node_pt(j)=old_node_pt;
1752  }
1753  }
1754 
1755 
1756  // Brick vertex node in centroid of tet. Only built for first element.
1757  // Enumerated "13" in initial sketch.
1758  {
1759  unsigned j=26;
1760 
1761  // Always copied
1762  el_pt->node_pt(j)=centroid_node_pt;
1763  }
1764 
1765 
1766  // Internal brick node -- always built
1767  {
1768  unsigned j=13;
1769 
1770  // Always new
1771  {
1772  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1773  Node_pt.push_back(new_node_pt);
1774  Vector<double> s(3);
1775  Vector<double> s_tet(3);
1776  Vector<double> x_tet(3);
1777  el_pt->local_coordinate_of_node(j,s);
1778  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1779  tet_el_pt->interpolated_x(s_tet,x_tet);
1780  new_node_pt->x(0)=x_tet[0];
1781  new_node_pt->x(1)=x_tet[1];
1782  new_node_pt->x(2)=x_tet[2];
1783  }
1784  }
1785 
1786 
1787  // Brick edge node between brick nodes 0 and 2
1788  {
1789  unsigned j=1;
1790 
1791  // Need new node?
1792  Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
1793  old_node_pt=brick_edge_node_pt[edge];
1794  if (old_node_pt==0)
1795  {
1796  Node* new_node_pt=0;
1797  if (edge.is_boundary_edge())
1798  {
1799  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1800  }
1801  else
1802  {
1803  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1804  }
1805  brick_edge_node_pt[edge]=new_node_pt;
1806  Node_pt.push_back(new_node_pt);
1807  Vector<double> s(3);
1808  Vector<double> s_tet(3);
1809  Vector<double> x_tet(3);
1810  el_pt->local_coordinate_of_node(j,s);
1811  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1812  tet_el_pt->interpolated_x(s_tet,x_tet);
1813  new_node_pt->x(0)=x_tet[0];
1814  new_node_pt->x(1)=x_tet[1];
1815  new_node_pt->x(2)=x_tet[2];
1816  }
1817  // Node already exists
1818  else
1819  {
1820  el_pt->node_pt(j)=old_node_pt;
1821  }
1822  }
1823 
1824 
1825  // Brick edge node between brick nodes 0 and 6
1826  {
1827  unsigned j=3;
1828 
1829  // Need new node?
1830  Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
1831  old_node_pt=brick_edge_node_pt[edge];
1832  if (old_node_pt==0)
1833  {
1834  Node* new_node_pt=0;
1835  if (edge.is_boundary_edge())
1836  {
1837  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1838  }
1839  else
1840  {
1841  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1842  }
1843  brick_edge_node_pt[edge]=new_node_pt;
1844  Node_pt.push_back(new_node_pt);
1845  Vector<double> s(3);
1846  Vector<double> s_tet(3);
1847  Vector<double> x_tet(3);
1848  el_pt->local_coordinate_of_node(j,s);
1849  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1850  tet_el_pt->interpolated_x(s_tet,x_tet);
1851  new_node_pt->x(0)=x_tet[0];
1852  new_node_pt->x(1)=x_tet[1];
1853  new_node_pt->x(2)=x_tet[2];
1854  }
1855  // Node already exists
1856  else
1857  {
1858  el_pt->node_pt(j)=old_node_pt;
1859  }
1860  }
1861 
1862 
1863  // Brick edge node between brick nodes 2 and 8
1864  {
1865  unsigned j=5;
1866 
1867  // Need new node?
1868  Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
1869  old_node_pt=brick_edge_node_pt[edge];
1870  if (old_node_pt==0)
1871  {
1872  Node* new_node_pt=0;
1873  if (edge.is_boundary_edge())
1874  {
1875  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1876  }
1877  else
1878  {
1879  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1880  }
1881  brick_edge_node_pt[edge]=new_node_pt;
1882  Node_pt.push_back(new_node_pt);
1883  Vector<double> s(3);
1884  Vector<double> s_tet(3);
1885  Vector<double> x_tet(3);
1886  el_pt->local_coordinate_of_node(j,s);
1887  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1888  tet_el_pt->interpolated_x(s_tet,x_tet);
1889  new_node_pt->x(0)=x_tet[0];
1890  new_node_pt->x(1)=x_tet[1];
1891  new_node_pt->x(2)=x_tet[2];
1892  }
1893  // Node already exists
1894  else
1895  {
1896  el_pt->node_pt(j)=old_node_pt;
1897  }
1898  }
1899 
1900  // Brick edge node between brick nodes 6 and 8
1901  {
1902  unsigned j=7;
1903 
1904  // Need new node?
1905  Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
1906  old_node_pt=brick_edge_node_pt[edge];
1907  if (old_node_pt==0)
1908  {
1909  Node* new_node_pt=0;
1910  if (edge.is_boundary_edge())
1911  {
1912  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1913  }
1914  else
1915  {
1916  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1917  }
1918  brick_edge_node_pt[edge]=new_node_pt;
1919  Node_pt.push_back(new_node_pt);
1920  Vector<double> s(3);
1921  Vector<double> s_tet(3);
1922  Vector<double> x_tet(3);
1923  el_pt->local_coordinate_of_node(j,s);
1924  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1925  tet_el_pt->interpolated_x(s_tet,x_tet);
1926  new_node_pt->x(0)=x_tet[0];
1927  new_node_pt->x(1)=x_tet[1];
1928  new_node_pt->x(2)=x_tet[2];
1929  }
1930  // Node already exists
1931  else
1932  {
1933  el_pt->node_pt(j)=old_node_pt;
1934  }
1935  }
1936 
1937  // Brick edge node between brick nodes 18 and 20
1938  {
1939  unsigned j=19;
1940 
1941  // Need new node?
1942  Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
1943  old_node_pt=brick_edge_node_pt[edge];
1944  if (old_node_pt==0)
1945  {
1946  Node* new_node_pt=0;
1947  if (edge.is_boundary_edge())
1948  {
1949  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1950  }
1951  else
1952  {
1953  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1954  }
1955  brick_edge_node_pt[edge]=new_node_pt;
1956  Node_pt.push_back(new_node_pt);
1957  Vector<double> s(3);
1958  Vector<double> s_tet(3);
1959  Vector<double> x_tet(3);
1960  el_pt->local_coordinate_of_node(j,s);
1961  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1962  tet_el_pt->interpolated_x(s_tet,x_tet);
1963  new_node_pt->x(0)=x_tet[0];
1964  new_node_pt->x(1)=x_tet[1];
1965  new_node_pt->x(2)=x_tet[2];
1966  }
1967  // Node already exists
1968  else
1969  {
1970  el_pt->node_pt(j)=old_node_pt;
1971  }
1972  }
1973 
1974 
1975  // Brick edge node between brick nodes 18 and 24
1976  {
1977  unsigned j=21;
1978 
1979  // Need new node?
1980  Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
1981  old_node_pt=brick_edge_node_pt[edge];
1982  if (old_node_pt==0)
1983  {
1984  Node* new_node_pt=0;
1985  if (edge.is_boundary_edge())
1986  {
1987  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1988  }
1989  else
1990  {
1991  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1992  }
1993  brick_edge_node_pt[edge]=new_node_pt;
1994  Node_pt.push_back(new_node_pt);
1995  Vector<double> s(3);
1996  Vector<double> s_tet(3);
1997  Vector<double> x_tet(3);
1998  el_pt->local_coordinate_of_node(j,s);
1999  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2000  tet_el_pt->interpolated_x(s_tet,x_tet);
2001  new_node_pt->x(0)=x_tet[0];
2002  new_node_pt->x(1)=x_tet[1];
2003  new_node_pt->x(2)=x_tet[2];
2004  }
2005  // Node already exists
2006  else
2007  {
2008  el_pt->node_pt(j)=old_node_pt;
2009  }
2010  }
2011 
2012  // Brick edge node between brick nodes 20 and 26
2013  {
2014  unsigned j=23;
2015 
2016  // Need new node?
2017  Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
2018  old_node_pt=brick_edge_node_pt[edge];
2019  if (old_node_pt==0)
2020  {
2021  Node* new_node_pt=0;
2022  if (edge.is_boundary_edge())
2023  {
2024  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2025  }
2026  else
2027  {
2028  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2029  }
2030  brick_edge_node_pt[edge]=new_node_pt;
2031  Node_pt.push_back(new_node_pt);
2032  Vector<double> s(3);
2033  Vector<double> s_tet(3);
2034  Vector<double> x_tet(3);
2035  el_pt->local_coordinate_of_node(j,s);
2036  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2037  tet_el_pt->interpolated_x(s_tet,x_tet);
2038  new_node_pt->x(0)=x_tet[0];
2039  new_node_pt->x(1)=x_tet[1];
2040  new_node_pt->x(2)=x_tet[2];
2041  }
2042  // Node already exists
2043  else
2044  {
2045  el_pt->node_pt(j)=old_node_pt;
2046  }
2047  }
2048 
2049 
2050  // Brick edge node between brick nodes 24 and 26
2051  {
2052  unsigned j=25;
2053 
2054  // Need new node?
2055  Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
2056  old_node_pt=brick_edge_node_pt[edge];
2057  if (old_node_pt==0)
2058  {
2059  Node* new_node_pt=0;
2060  if (edge.is_boundary_edge())
2061  {
2062  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2063  }
2064  else
2065  {
2066  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2067  }
2068  brick_edge_node_pt[edge]=new_node_pt;
2069  Node_pt.push_back(new_node_pt);
2070  Vector<double> s(3);
2071  Vector<double> s_tet(3);
2072  Vector<double> x_tet(3);
2073  el_pt->local_coordinate_of_node(j,s);
2074  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2075  tet_el_pt->interpolated_x(s_tet,x_tet);
2076  new_node_pt->x(0)=x_tet[0];
2077  new_node_pt->x(1)=x_tet[1];
2078  new_node_pt->x(2)=x_tet[2];
2079  }
2080  // Node already exists
2081  else
2082  {
2083  el_pt->node_pt(j)=old_node_pt;
2084  }
2085  }
2086 
2087  // Brick edge node between brick nodes 0 and 18
2088  {
2089  unsigned j=9;
2090 
2091  // Need new node?
2092  Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
2093  old_node_pt=brick_edge_node_pt[edge];
2094  if (old_node_pt==0)
2095  {
2096  Node* new_node_pt=0;
2097  if (edge.is_boundary_edge())
2098  {
2099  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2100  }
2101  else
2102  {
2103  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2104  }
2105  brick_edge_node_pt[edge]=new_node_pt;
2106  Node_pt.push_back(new_node_pt);
2107  Vector<double> s(3);
2108  Vector<double> s_tet(3);
2109  Vector<double> x_tet(3);
2110  el_pt->local_coordinate_of_node(j,s);
2111  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2112  tet_el_pt->interpolated_x(s_tet,x_tet);
2113  new_node_pt->x(0)=x_tet[0];
2114  new_node_pt->x(1)=x_tet[1];
2115  new_node_pt->x(2)=x_tet[2];
2116  }
2117  // Node already exists
2118  else
2119  {
2120  el_pt->node_pt(j)=old_node_pt;
2121  }
2122  }
2123 
2124 
2125  // Brick edge node between brick nodes 2 and 20
2126  {
2127  unsigned j=11;
2128 
2129  // Need new node?
2130  Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
2131  old_node_pt=brick_edge_node_pt[edge];
2132  if (old_node_pt==0)
2133  {
2134  Node* new_node_pt=0;
2135  if (edge.is_boundary_edge())
2136  {
2137  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2138  }
2139  else
2140  {
2141  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2142  }
2143  brick_edge_node_pt[edge]=new_node_pt;
2144  Node_pt.push_back(new_node_pt);
2145  Vector<double> s(3);
2146  Vector<double> s_tet(3);
2147  Vector<double> x_tet(3);
2148  el_pt->local_coordinate_of_node(j,s);
2149  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2150  tet_el_pt->interpolated_x(s_tet,x_tet);
2151  new_node_pt->x(0)=x_tet[0];
2152  new_node_pt->x(1)=x_tet[1];
2153  new_node_pt->x(2)=x_tet[2];
2154  }
2155  // Node already exists
2156  else
2157  {
2158  el_pt->node_pt(j)=old_node_pt;
2159  }
2160  }
2161 
2162 
2163  // Brick edge node between brick nodes 6 and 24
2164  {
2165  unsigned j=15;
2166 
2167  // Need new node?
2168  Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
2169  old_node_pt=brick_edge_node_pt[edge];
2170  if (old_node_pt==0)
2171  {
2172  Node* new_node_pt=0;
2173  if (edge.is_boundary_edge())
2174  {
2175  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2176  }
2177  else
2178  {
2179  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2180  }
2181  brick_edge_node_pt[edge]=new_node_pt;
2182  Node_pt.push_back(new_node_pt);
2183  Vector<double> s(3);
2184  Vector<double> s_tet(3);
2185  Vector<double> x_tet(3);
2186  el_pt->local_coordinate_of_node(j,s);
2187  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2188  tet_el_pt->interpolated_x(s_tet,x_tet);
2189  new_node_pt->x(0)=x_tet[0];
2190  new_node_pt->x(1)=x_tet[1];
2191  new_node_pt->x(2)=x_tet[2];
2192  }
2193  // Node already exists
2194  else
2195  {
2196  el_pt->node_pt(j)=old_node_pt;
2197  }
2198  }
2199 
2200 
2201  // Brick edge node between brick nodes 8 and 26
2202  {
2203  unsigned j=17;
2204 
2205  // Need new node?
2206  Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
2207  old_node_pt=brick_edge_node_pt[edge];
2208  if (old_node_pt==0)
2209  {
2210  Node* new_node_pt=0;
2211  if (edge.is_boundary_edge())
2212  {
2213  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2214  }
2215  else
2216  {
2217  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2218  }
2219  brick_edge_node_pt[edge]=new_node_pt;
2220  Node_pt.push_back(new_node_pt);
2221  Vector<double> s(3);
2222  Vector<double> s_tet(3);
2223  Vector<double> x_tet(3);
2224  el_pt->local_coordinate_of_node(j,s);
2225  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2226  tet_el_pt->interpolated_x(s_tet,x_tet);
2227  new_node_pt->x(0)=x_tet[0];
2228  new_node_pt->x(1)=x_tet[1];
2229  new_node_pt->x(2)=x_tet[2];
2230  }
2231  // Node already exists
2232  else
2233  {
2234  el_pt->node_pt(j)=old_node_pt;
2235  }
2236  }
2237 
2238 
2239  // Mid brick-face node associated with face
2240  // spanned by mid-vertex nodes associated with tet edges 0 and 2
2241  {
2242  unsigned j=10;
2243 
2244  // Need new node?
2245  TFace face(tet_el_pt->node_pt( central_tet_vertex),
2246  tet_el_pt->node_pt(tet_edge_node[0][0]),
2247  tet_el_pt->node_pt(tet_edge_node[2][0]));
2248 
2249  old_node_pt=tet_face_node_pt[face];
2250  if (old_node_pt==0)
2251  {
2252  Node* new_node_pt=0;
2253  if (face.is_boundary_face())
2254  {
2255  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2256  }
2257  else
2258  {
2259  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2260  }
2261  tet_face_node_pt[face]=new_node_pt;
2262  Node_pt.push_back(new_node_pt);
2263  Vector<double> s(3);
2264  Vector<double> s_tet(3);
2265  Vector<double> x_tet(3);
2266  el_pt->local_coordinate_of_node(j,s);
2267  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2268  tet_el_pt->interpolated_x(s_tet,x_tet);
2269  new_node_pt->x(0)=x_tet[0];
2270  new_node_pt->x(1)=x_tet[1];
2271  new_node_pt->x(2)=x_tet[2];
2272  }
2273  // Node already exists
2274  else
2275  {
2276  el_pt->node_pt(j)=old_node_pt;
2277  }
2278  }
2279 
2280 
2281 
2282  // Mid brick-face node associated with face
2283  // spanned by mid-vertex nodes associated with tet edges 1 and 2
2284  {
2285  unsigned j=12;
2286 
2287  // Need new node?
2288  TFace face(tet_el_pt->node_pt(central_tet_vertex),
2289  tet_el_pt->node_pt(tet_edge_node[1][0]),
2290  tet_el_pt->node_pt(tet_edge_node[2][0]));
2291 
2292  old_node_pt=tet_face_node_pt[face];
2293  if (old_node_pt==0)
2294  {
2295  Node* new_node_pt=0;
2296  if (face.is_boundary_face())
2297  {
2298  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2299  }
2300  else
2301  {
2302  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2303  }
2304  tet_face_node_pt[face]=new_node_pt;
2305  Node_pt.push_back(new_node_pt);
2306  Vector<double> s(3);
2307  Vector<double> s_tet(3);
2308  Vector<double> x_tet(3);
2309  el_pt->local_coordinate_of_node(j,s);
2310  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2311  tet_el_pt->interpolated_x(s_tet,x_tet);
2312  new_node_pt->x(0)=x_tet[0];
2313  new_node_pt->x(1)=x_tet[1];
2314  new_node_pt->x(2)=x_tet[2];
2315  }
2316  // Node already exists
2317  else
2318  {
2319  el_pt->node_pt(j)=old_node_pt;
2320  }
2321  }
2322 
2323 
2324 
2325  // Mid brick-face node associated with face
2326  // spanned by mid-vertex nodes associated with tet edges 0 and 1
2327  {
2328  unsigned j=4;
2329 
2330  // Need new node?
2331  TFace face(tet_el_pt->node_pt(central_tet_vertex),
2332  tet_el_pt->node_pt(tet_edge_node[0][0]),
2333  tet_el_pt->node_pt(tet_edge_node[1][0]));
2334 
2335  old_node_pt=tet_face_node_pt[face];
2336  if (old_node_pt==0)
2337  {
2338  Node* new_node_pt=0;
2339  if (face.is_boundary_face())
2340  {
2341  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2342  }
2343  else
2344  {
2345  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2346  }
2347  tet_face_node_pt[face]=new_node_pt;
2348  Node_pt.push_back(new_node_pt);
2349  Vector<double> s(3);
2350  Vector<double> s_tet(3);
2351  Vector<double> x_tet(3);
2352  el_pt->local_coordinate_of_node(j,s);
2353  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2354  tet_el_pt->interpolated_x(s_tet,x_tet);
2355  new_node_pt->x(0)=x_tet[0];
2356  new_node_pt->x(1)=x_tet[1];
2357  new_node_pt->x(2)=x_tet[2];
2358  }
2359  // Node already exists
2360  else
2361  {
2362  el_pt->node_pt(j)=old_node_pt;
2363  }
2364  }
2365 
2366 
2367  // Top mid brick-face node -- only built by this element
2368  {
2369  unsigned j=22;
2370  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2371  Node_pt.push_back(new_node_pt);
2372  Vector<double> s(3);
2373  Vector<double> s_tet(3);
2374  Vector<double> x_tet(3);
2375  el_pt->local_coordinate_of_node(j,s);
2376  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2377  top_mid_face_node1_pt=new_node_pt;
2378  tet_el_pt->interpolated_x(s_tet,x_tet);
2379  new_node_pt->x(0)=x_tet[0];
2380  new_node_pt->x(1)=x_tet[1];
2381  new_node_pt->x(2)=x_tet[2];
2382  }
2383 
2384 
2385 
2386  // Right mid brick-face node -- only built by this element
2387  {
2388  unsigned j=14;
2389  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2390  Node_pt.push_back(new_node_pt);
2391  Vector<double> s(3);
2392  Vector<double> s_tet(3);
2393  Vector<double> x_tet(3);
2394  el_pt->local_coordinate_of_node(j,s);
2395  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2396  right_mid_face_node1_pt=new_node_pt;
2397  tet_el_pt->interpolated_x(s_tet,x_tet);
2398  new_node_pt->x(0)=x_tet[0];
2399  new_node_pt->x(1)=x_tet[1];
2400  new_node_pt->x(2)=x_tet[2];
2401  }
2402 
2403 
2404  // Back mid brick-face node copy from previous element
2405  {
2406  unsigned j=16;
2407 
2408  // Always copied
2409  el_pt->node_pt(j)=right_mid_face_node0_pt;
2410  }
2411  }
2412 
2413 
2414  // Third brick element is centred at node 3 of tet:
2415  //-------------------------------------------------
2416  {
2417  // Assign coordinates of dummy element
2418  for (unsigned j=0;j<8;j++)
2419  {
2420  Node* nod_pt=dummy_q_el_pt[2]->node_pt(j);
2421  Vector<double> s_tet(3);
2422  Vector<double> x_tet(3);
2423  switch (j)
2424  {
2425  case 0:
2426  tet_el_pt->local_coordinate_of_node(3,s_tet);
2427  nod_pt->set_value(0,s_tet[0]);
2428  nod_pt->set_value(1,s_tet[1]);
2429  nod_pt->set_value(2,s_tet[2]);
2430  tet_el_pt->interpolated_x(s_tet,x_tet);
2431  nod_pt->x(0)=x_tet[0];
2432  nod_pt->x(1)=x_tet[1];
2433  nod_pt->x(2)=x_tet[2];
2434  break;
2435  case 1:
2436  tet_el_pt->local_coordinate_of_node(6,s_tet);
2437  nod_pt->set_value(0,s_tet[0]);
2438  nod_pt->set_value(1,s_tet[1]);
2439  nod_pt->set_value(2,s_tet[2]);
2440  tet_el_pt->interpolated_x(s_tet,x_tet);
2441  nod_pt->x(0)=x_tet[0];
2442  nod_pt->x(1)=x_tet[1];
2443  nod_pt->x(2)=x_tet[2];
2444  break;
2445  case 2:
2446  tet_el_pt->local_coordinate_of_node(9,s_tet);
2447  nod_pt->set_value(0,s_tet[0]);
2448  nod_pt->set_value(1,s_tet[1]);
2449  nod_pt->set_value(2,s_tet[2]);
2450  tet_el_pt->interpolated_x(s_tet,x_tet);
2451  nod_pt->x(0)=x_tet[0];
2452  nod_pt->x(1)=x_tet[1];
2453  nod_pt->x(2)=x_tet[2];
2454  break;
2455  case 3:
2456  // label 13 in initial sketch: Mid face node on face
2457  // spanned by tet nodes 0,1,3
2458  s_tet[0]=1.0/3.0;
2459  s_tet[1]=1.0/3.0;
2460  s_tet[2]=0.0;
2461  nod_pt->set_value(0,s_tet[0]);
2462  nod_pt->set_value(1,s_tet[1]);
2463  nod_pt->set_value(2,s_tet[2]);
2464  tet_el_pt->interpolated_x(s_tet,x_tet);
2465  nod_pt->x(0)=x_tet[0];
2466  nod_pt->x(1)=x_tet[1];
2467  nod_pt->x(2)=x_tet[2];
2468  break;
2469  case 4:
2470  tet_el_pt->local_coordinate_of_node(8,s_tet);
2471  nod_pt->set_value(0,s_tet[0]);
2472  nod_pt->set_value(1,s_tet[1]);
2473  nod_pt->set_value(2,s_tet[2]);
2474  tet_el_pt->interpolated_x(s_tet,x_tet);
2475  nod_pt->x(0)=x_tet[0];
2476  nod_pt->x(1)=x_tet[1];
2477  nod_pt->x(2)=x_tet[2];
2478  break;
2479  case 5:
2480  // label 12 in initial sketch: Mid face node on face
2481  // spanned by tet nodes 0,2,3
2482  s_tet[0]=1.0/3.0;
2483  s_tet[1]=0.0;
2484  s_tet[2]=1.0/3.0;
2485  nod_pt->set_value(0,s_tet[0]);
2486  nod_pt->set_value(1,s_tet[1]);
2487  nod_pt->set_value(2,s_tet[2]);
2488  tet_el_pt->interpolated_x(s_tet,x_tet);
2489  nod_pt->x(0)=x_tet[0];
2490  nod_pt->x(1)=x_tet[1];
2491  nod_pt->x(2)=x_tet[2];
2492  break;
2493  case 6:
2494  // label 10 in initial sketch: Mid face node on face
2495  // spanned by tet nodes 1,2,3
2496  s_tet[0]=0.0;
2497  s_tet[1]=1.0/3.0;
2498  s_tet[2]=1.0/3.0;
2499  nod_pt->set_value(0,s_tet[0]);
2500  nod_pt->set_value(1,s_tet[1]);
2501  nod_pt->set_value(2,s_tet[2]);
2502  tet_el_pt->interpolated_x(s_tet,x_tet);
2503  nod_pt->x(0)=x_tet[0];
2504  nod_pt->x(1)=x_tet[1];
2505  nod_pt->x(2)=x_tet[2];
2506  break;
2507  case 7:
2508  // label 14 in initial sketch: Centroid
2509  s_tet[0]=0.25;
2510  s_tet[1]=0.25;
2511  s_tet[2]=0.25;
2512  nod_pt->set_value(0,s_tet[0]);
2513  nod_pt->set_value(1,s_tet[1]);
2514  nod_pt->set_value(2,s_tet[2]);
2515  tet_el_pt->interpolated_x(s_tet,x_tet);
2516  nod_pt->x(0)=x_tet[0];
2517  nod_pt->x(1)=x_tet[1];
2518  nod_pt->x(2)=x_tet[2];
2519  break;
2520  }
2521  }
2522 
2523 
2524  // Create actual second brick element
2525  FiniteElement* el_pt=new ELEMENT;
2526  brick_el2_pt=el_pt;
2527  Element_pt.push_back(el_pt);
2528 
2529  TFace face0(tet_el_pt->node_pt(3),
2530  tet_el_pt->node_pt(0),
2531  tet_el_pt->node_pt(2));
2532 
2533  TFace face1(tet_el_pt->node_pt(3),
2534  tet_el_pt->node_pt(1),
2535  tet_el_pt->node_pt(2));
2536 
2537  TFace face2(tet_el_pt->node_pt(3),
2538  tet_el_pt->node_pt(1),
2539  tet_el_pt->node_pt(0));
2540 
2541  // Tet vertex nodes along edges emanating from node 0 in brick
2542  Vector<Vector<unsigned> > tet_edge_node(3);
2543  tet_edge_node[0].resize(2);
2544  tet_edge_node[0][0]=6;
2545  tet_edge_node[0][1]=0;
2546  tet_edge_node[1].resize(2);
2547  tet_edge_node[1][0]=9;
2548  tet_edge_node[1][1]=1;
2549  tet_edge_node[2].resize(2);
2550  tet_edge_node[2][0]=8;
2551  tet_edge_node[2][1]=2;
2552 
2553  // Node number of tet vertex that node 0 in brick is centred on
2554  unsigned central_tet_vertex=3;
2555 
2556  Node* tet_node_pt=0;
2557  Node* old_node_pt=0;
2558 
2559  // Corner node
2560  {
2561  unsigned j=0;
2562 
2563  // Need new node?
2564  tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
2565  old_node_pt=tet_node_node_pt[tet_node_pt];
2566  if (old_node_pt==0)
2567  {
2568  Node* new_node_pt=0;
2569  if (tet_node_pt->is_on_boundary())
2570  {
2571  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2572  }
2573  else
2574  {
2575  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2576  }
2577  tet_node_node_pt[tet_node_pt]=new_node_pt;
2578  Node_pt.push_back(new_node_pt);
2579  Vector<double> s(3);
2580  Vector<double> s_tet(3);
2581  Vector<double> x_tet(3);
2582  el_pt->local_coordinate_of_node(j,s);
2583  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2584  tet_el_pt->interpolated_x(s_tet,x_tet);
2585  new_node_pt->x(0)=x_tet[0];
2586  new_node_pt->x(1)=x_tet[1];
2587  new_node_pt->x(2)=x_tet[2];
2588  }
2589  // Node already exists
2590  else
2591  {
2592  el_pt->node_pt(j)=old_node_pt;
2593  }
2594  }
2595 
2596 
2597  // Brick vertex node coindides with mid-edge node on tet edge 0
2598  {
2599  unsigned j=2;
2600 
2601  // Need new node?
2602  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
2603  old_node_pt=tet_node_node_pt[tet_node_pt];
2604  if (old_node_pt==0)
2605  {
2606  Node* new_node_pt=0;
2607  if (tet_node_pt->is_on_boundary())
2608  {
2609  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2610  }
2611  else
2612  {
2613  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2614  }
2615  tet_node_node_pt[tet_node_pt]=new_node_pt;
2616  Node_pt.push_back(new_node_pt);
2617  Vector<double> s(3);
2618  Vector<double> s_tet(3);
2619  Vector<double> x_tet(3);
2620  el_pt->local_coordinate_of_node(j,s);
2621  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2622  tet_el_pt->interpolated_x(s_tet,x_tet);
2623  new_node_pt->x(0)=x_tet[0];
2624  new_node_pt->x(1)=x_tet[1];
2625  new_node_pt->x(2)=x_tet[2];
2626  }
2627  // Node already exists
2628  else
2629  {
2630  el_pt->node_pt(j)=old_node_pt;
2631  }
2632  }
2633 
2634 
2635  // Brick vertex node coindides with mid vertex node of tet edge 1
2636  {
2637  unsigned j=6;
2638 
2639  // Need new node?
2640  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
2641  old_node_pt=tet_node_node_pt[tet_node_pt];
2642  if (old_node_pt==0)
2643  {
2644  Node* new_node_pt=0;
2645  if (tet_node_pt->is_on_boundary())
2646  {
2647  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2648  }
2649  else
2650  {
2651  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2652  }
2653  tet_node_node_pt[tet_node_pt]=new_node_pt;
2654  Node_pt.push_back(new_node_pt);
2655  Vector<double> s(3);
2656  Vector<double> s_tet(3);
2657  Vector<double> x_tet(3);
2658  el_pt->local_coordinate_of_node(j,s);
2659  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2660  tet_el_pt->interpolated_x(s_tet,x_tet);
2661  new_node_pt->x(0)=x_tet[0];
2662  new_node_pt->x(1)=x_tet[1];
2663  new_node_pt->x(2)=x_tet[2];
2664  }
2665  // Node already exists
2666  else
2667  {
2668  el_pt->node_pt(j)=old_node_pt;
2669  }
2670  }
2671 
2672 
2673  // Brick vertex node coindides with mid-vertex node of tet edge 2
2674  {
2675  unsigned j=18;
2676 
2677  // Need new node?
2678  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
2679  old_node_pt=tet_node_node_pt[tet_node_pt];
2680  if (old_node_pt==0)
2681  {
2682  Node* new_node_pt=0;
2683  if (tet_node_pt->is_on_boundary())
2684  {
2685  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2686  }
2687  else
2688  {
2689  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2690  }
2691  tet_node_node_pt[tet_node_pt]=new_node_pt;
2692  Node_pt.push_back(new_node_pt);
2693  Vector<double> s(3);
2694  Vector<double> s_tet(3);
2695  Vector<double> x_tet(3);
2696  el_pt->local_coordinate_of_node(j,s);
2697  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2698  tet_el_pt->interpolated_x(s_tet,x_tet);
2699  new_node_pt->x(0)=x_tet[0];
2700  new_node_pt->x(1)=x_tet[1];
2701  new_node_pt->x(2)=x_tet[2];
2702  }
2703  // Node already exists
2704  else
2705  {
2706  el_pt->node_pt(j)=old_node_pt;
2707  }
2708  }
2709 
2710 
2711 
2712  // Brick vertex node in the middle of tet face0
2713  {
2714  unsigned j=20;
2715 
2716  // Need new node?
2717  old_node_pt=tet_face_node_pt[face0];
2718  if (old_node_pt==0)
2719  {
2720  Node* new_node_pt=0;
2721  if (face0.is_boundary_face())
2722  {
2723  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2724  }
2725  else
2726  {
2727  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2728  }
2729  tet_face_node_pt[face0]=new_node_pt;
2730  Node_pt.push_back(new_node_pt);
2731  Vector<double> s(3);
2732  Vector<double> s_tet(3);
2733  Vector<double> x_tet(3);
2734  el_pt->local_coordinate_of_node(j,s);
2735  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2736  tet_el_pt->interpolated_x(s_tet,x_tet);
2737  new_node_pt->x(0)=x_tet[0];
2738  new_node_pt->x(1)=x_tet[1];
2739  new_node_pt->x(2)=x_tet[2];
2740  }
2741  // Node already exists
2742  else
2743  {
2744  el_pt->node_pt(j)=old_node_pt;
2745  }
2746  }
2747 
2748  // Brick vertex node in the middle of tet face1
2749  {
2750  unsigned j=24;
2751 
2752  // Need new node?
2753  old_node_pt=tet_face_node_pt[face1];
2754  if (old_node_pt==0)
2755  {
2756  Node* new_node_pt=0;
2757  if (face1.is_boundary_face())
2758  {
2759  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2760  }
2761  else
2762  {
2763  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2764  }
2765  tet_face_node_pt[face1]=new_node_pt;
2766  Node_pt.push_back(new_node_pt);
2767  Vector<double> s(3);
2768  Vector<double> s_tet(3);
2769  Vector<double> x_tet(3);
2770  el_pt->local_coordinate_of_node(j,s);
2771  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2772  tet_el_pt->interpolated_x(s_tet,x_tet);
2773  new_node_pt->x(0)=x_tet[0];
2774  new_node_pt->x(1)=x_tet[1];
2775  new_node_pt->x(2)=x_tet[2];
2776  }
2777  // Node already exists
2778  else
2779  {
2780  el_pt->node_pt(j)=old_node_pt;
2781  }
2782  }
2783 
2784  // Brick vertex node in the middle of tet face2
2785  {
2786  unsigned j=8;
2787 
2788  // Need new node?
2789  old_node_pt=tet_face_node_pt[face2];
2790  if (old_node_pt==0)
2791  {
2792  Node* new_node_pt=0;
2793  if (face2.is_boundary_face())
2794  {
2795  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2796  }
2797  else
2798  {
2799  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2800  }
2801  tet_face_node_pt[face2]=new_node_pt;
2802  Node_pt.push_back(new_node_pt);
2803  Vector<double> s(3);
2804  Vector<double> s_tet(3);
2805  Vector<double> x_tet(3);
2806  el_pt->local_coordinate_of_node(j,s);
2807  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2808  tet_el_pt->interpolated_x(s_tet,x_tet);
2809  new_node_pt->x(0)=x_tet[0];
2810  new_node_pt->x(1)=x_tet[1];
2811  new_node_pt->x(2)=x_tet[2];
2812  }
2813  // Node already exists
2814  else
2815  {
2816  el_pt->node_pt(j)=old_node_pt;
2817  }
2818  }
2819 
2820 
2821  // Brick vertex node in centroid of tet. Only built for first element.
2822  // Enumerated "13" in initial sketch.
2823  {
2824  unsigned j=26;
2825 
2826  // Always copied
2827  el_pt->node_pt(j)=centroid_node_pt;
2828  }
2829 
2830 
2831  // Internal brick node -- always built
2832  {
2833  unsigned j=13;
2834 
2835  // Always new
2836  {
2837  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2838  Node_pt.push_back(new_node_pt);
2839  Vector<double> s(3);
2840  Vector<double> s_tet(3);
2841  Vector<double> x_tet(3);
2842  el_pt->local_coordinate_of_node(j,s);
2843  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2844  tet_el_pt->interpolated_x(s_tet,x_tet);
2845  new_node_pt->x(0)=x_tet[0];
2846  new_node_pt->x(1)=x_tet[1];
2847  new_node_pt->x(2)=x_tet[2];
2848  }
2849  }
2850 
2851  // Brick edge node between brick nodes 0 and 2
2852  {
2853  unsigned j=1;
2854 
2855  // Need new node?
2856  Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
2857  old_node_pt=brick_edge_node_pt[edge];
2858  if (old_node_pt==0)
2859  {
2860  Node* new_node_pt=0;
2861  if (edge.is_boundary_edge())
2862  {
2863  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2864  }
2865  else
2866  {
2867  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2868  }
2869  brick_edge_node_pt[edge]=new_node_pt;
2870  Node_pt.push_back(new_node_pt);
2871  Vector<double> s(3);
2872  Vector<double> s_tet(3);
2873  Vector<double> x_tet(3);
2874  el_pt->local_coordinate_of_node(j,s);
2875  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2876  tet_el_pt->interpolated_x(s_tet,x_tet);
2877  new_node_pt->x(0)=x_tet[0];
2878  new_node_pt->x(1)=x_tet[1];
2879  new_node_pt->x(2)=x_tet[2];
2880  }
2881  // Node already exists
2882  else
2883  {
2884  el_pt->node_pt(j)=old_node_pt;
2885  }
2886  }
2887 
2888 
2889  // Brick edge node between brick nodes 0 and 6
2890  {
2891  unsigned j=3;
2892 
2893  // Need new node?
2894  Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
2895  old_node_pt=brick_edge_node_pt[edge];
2896  if (old_node_pt==0)
2897  {
2898  Node* new_node_pt=0;
2899  if (edge.is_boundary_edge())
2900  {
2901  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2902  }
2903  else
2904  {
2905  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2906  }
2907  brick_edge_node_pt[edge]=new_node_pt;
2908  Node_pt.push_back(new_node_pt);
2909  Vector<double> s(3);
2910  Vector<double> s_tet(3);
2911  Vector<double> x_tet(3);
2912  el_pt->local_coordinate_of_node(j,s);
2913  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2914  tet_el_pt->interpolated_x(s_tet,x_tet);
2915  new_node_pt->x(0)=x_tet[0];
2916  new_node_pt->x(1)=x_tet[1];
2917  new_node_pt->x(2)=x_tet[2];
2918  }
2919  // Node already exists
2920  else
2921  {
2922  el_pt->node_pt(j)=old_node_pt;
2923  }
2924  }
2925 
2926  // Brick edge node between brick nodes 2 and 8
2927  {
2928  unsigned j=5;
2929 
2930  // Need new node?
2931  Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
2932  old_node_pt=brick_edge_node_pt[edge];
2933  if (old_node_pt==0)
2934  {
2935  Node* new_node_pt=0;
2936  if (edge.is_boundary_edge())
2937  {
2938  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2939  }
2940  else
2941  {
2942  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2943  }
2944  brick_edge_node_pt[edge]=new_node_pt;
2945  Node_pt.push_back(new_node_pt);
2946  Vector<double> s(3);
2947  Vector<double> s_tet(3);
2948  Vector<double> x_tet(3);
2949  el_pt->local_coordinate_of_node(j,s);
2950  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2951  tet_el_pt->interpolated_x(s_tet,x_tet);
2952  new_node_pt->x(0)=x_tet[0];
2953  new_node_pt->x(1)=x_tet[1];
2954  new_node_pt->x(2)=x_tet[2];
2955  }
2956  // Node already exists
2957  else
2958  {
2959  el_pt->node_pt(j)=old_node_pt;
2960  }
2961  }
2962 
2963  // Brick edge node between brick nodes 6 and 8
2964  {
2965  unsigned j=7;
2966 
2967  // Need new node?
2968  Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
2969  old_node_pt=brick_edge_node_pt[edge];
2970  if (old_node_pt==0)
2971  {
2972  Node* new_node_pt=0;
2973  if (edge.is_boundary_edge())
2974  {
2975  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2976  }
2977  else
2978  {
2979  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2980  }
2981  brick_edge_node_pt[edge]=new_node_pt;
2982  Node_pt.push_back(new_node_pt);
2983  Vector<double> s(3);
2984  Vector<double> s_tet(3);
2985  Vector<double> x_tet(3);
2986  el_pt->local_coordinate_of_node(j,s);
2987  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2988  tet_el_pt->interpolated_x(s_tet,x_tet);
2989  new_node_pt->x(0)=x_tet[0];
2990  new_node_pt->x(1)=x_tet[1];
2991  new_node_pt->x(2)=x_tet[2];
2992  }
2993  // Node already exists
2994  else
2995  {
2996  el_pt->node_pt(j)=old_node_pt;
2997  }
2998  }
2999 
3000  // Brick edge node between brick nodes 18 and 20
3001  {
3002  unsigned j=19;
3003 
3004  // Need new node?
3005  Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
3006  old_node_pt=brick_edge_node_pt[edge];
3007  if (old_node_pt==0)
3008  {
3009  Node* new_node_pt=0;
3010  if (edge.is_boundary_edge())
3011  {
3012  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3013  }
3014  else
3015  {
3016  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3017  }
3018  brick_edge_node_pt[edge]=new_node_pt;
3019  Node_pt.push_back(new_node_pt);
3020  Vector<double> s(3);
3021  Vector<double> s_tet(3);
3022  Vector<double> x_tet(3);
3023  el_pt->local_coordinate_of_node(j,s);
3024  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3025  tet_el_pt->interpolated_x(s_tet,x_tet);
3026  new_node_pt->x(0)=x_tet[0];
3027  new_node_pt->x(1)=x_tet[1];
3028  new_node_pt->x(2)=x_tet[2];
3029  }
3030  // Node already exists
3031  else
3032  {
3033  el_pt->node_pt(j)=old_node_pt;
3034  }
3035  }
3036 
3037 
3038  // Brick edge node between brick nodes 18 and 24
3039  {
3040  unsigned j=21;
3041 
3042  // Need new node?
3043  Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
3044  old_node_pt=brick_edge_node_pt[edge];
3045  if (old_node_pt==0)
3046  {
3047  Node* new_node_pt=0;
3048  if (edge.is_boundary_edge())
3049  {
3050  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3051  }
3052  else
3053  {
3054  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3055  }
3056  brick_edge_node_pt[edge]=new_node_pt;
3057  Node_pt.push_back(new_node_pt);
3058  Vector<double> s(3);
3059  Vector<double> s_tet(3);
3060  Vector<double> x_tet(3);
3061  el_pt->local_coordinate_of_node(j,s);
3062  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3063  tet_el_pt->interpolated_x(s_tet,x_tet);
3064  new_node_pt->x(0)=x_tet[0];
3065  new_node_pt->x(1)=x_tet[1];
3066  new_node_pt->x(2)=x_tet[2];
3067  }
3068  // Node already exists
3069  else
3070  {
3071  el_pt->node_pt(j)=old_node_pt;
3072  }
3073  }
3074 
3075  // Brick edge node between brick nodes 20 and 26
3076  {
3077  unsigned j=23;
3078 
3079  // Need new node?
3080  Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
3081  old_node_pt=brick_edge_node_pt[edge];
3082  if (old_node_pt==0)
3083  {
3084  Node* new_node_pt=0;
3085  if (edge.is_boundary_edge())
3086  {
3087  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3088  }
3089  else
3090  {
3091  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3092  }
3093  brick_edge_node_pt[edge]=new_node_pt;
3094  Node_pt.push_back(new_node_pt);
3095  Vector<double> s(3);
3096  Vector<double> s_tet(3);
3097  Vector<double> x_tet(3);
3098  el_pt->local_coordinate_of_node(j,s);
3099  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3100  tet_el_pt->interpolated_x(s_tet,x_tet);
3101  new_node_pt->x(0)=x_tet[0];
3102  new_node_pt->x(1)=x_tet[1];
3103  new_node_pt->x(2)=x_tet[2];
3104  }
3105  // Node already exists
3106  else
3107  {
3108  el_pt->node_pt(j)=old_node_pt;
3109  }
3110  }
3111 
3112 
3113  // Brick edge node between brick nodes 24 and 26
3114  {
3115  unsigned j=25;
3116 
3117  // Need new node?
3118  Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
3119  old_node_pt=brick_edge_node_pt[edge];
3120  if (old_node_pt==0)
3121  {
3122  Node* new_node_pt=0;
3123  if (edge.is_boundary_edge())
3124  {
3125  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3126  }
3127  else
3128  {
3129  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3130  }
3131  brick_edge_node_pt[edge]=new_node_pt;
3132  Node_pt.push_back(new_node_pt);
3133  Vector<double> s(3);
3134  Vector<double> s_tet(3);
3135  Vector<double> x_tet(3);
3136  el_pt->local_coordinate_of_node(j,s);
3137  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3138  tet_el_pt->interpolated_x(s_tet,x_tet);
3139  new_node_pt->x(0)=x_tet[0];
3140  new_node_pt->x(1)=x_tet[1];
3141  new_node_pt->x(2)=x_tet[2];
3142  }
3143  // Node already exists
3144  else
3145  {
3146  el_pt->node_pt(j)=old_node_pt;
3147  }
3148  }
3149 
3150  // Brick edge node between brick nodes 0 and 18
3151  {
3152  unsigned j=9;
3153 
3154  // Need new node?
3155  Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
3156  old_node_pt=brick_edge_node_pt[edge];
3157  if (old_node_pt==0)
3158  {
3159  Node* new_node_pt=0;
3160  if (edge.is_boundary_edge())
3161  {
3162  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3163  }
3164  else
3165  {
3166  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3167  }
3168  brick_edge_node_pt[edge]=new_node_pt;
3169  Node_pt.push_back(new_node_pt);
3170  Vector<double> s(3);
3171  Vector<double> s_tet(3);
3172  Vector<double> x_tet(3);
3173  el_pt->local_coordinate_of_node(j,s);
3174  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3175  tet_el_pt->interpolated_x(s_tet,x_tet);
3176  new_node_pt->x(0)=x_tet[0];
3177  new_node_pt->x(1)=x_tet[1];
3178  new_node_pt->x(2)=x_tet[2];
3179  }
3180  // Node already exists
3181  else
3182  {
3183  el_pt->node_pt(j)=old_node_pt;
3184  }
3185  }
3186 
3187 
3188  // Brick edge node between brick nodes 2 and 20
3189  {
3190  unsigned j=11;
3191 
3192  // Need new node?
3193  Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
3194  old_node_pt=brick_edge_node_pt[edge];
3195  if (old_node_pt==0)
3196  {
3197  Node* new_node_pt=0;
3198  if (edge.is_boundary_edge())
3199  {
3200  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3201  }
3202  else
3203  {
3204  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3205  }
3206  brick_edge_node_pt[edge]=new_node_pt;
3207  Node_pt.push_back(new_node_pt);
3208  Vector<double> s(3);
3209  Vector<double> s_tet(3);
3210  Vector<double> x_tet(3);
3211  el_pt->local_coordinate_of_node(j,s);
3212  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3213  tet_el_pt->interpolated_x(s_tet,x_tet);
3214  new_node_pt->x(0)=x_tet[0];
3215  new_node_pt->x(1)=x_tet[1];
3216  new_node_pt->x(2)=x_tet[2];
3217  }
3218  // Node already exists
3219  else
3220  {
3221  el_pt->node_pt(j)=old_node_pt;
3222  }
3223  }
3224 
3225 
3226  // Brick edge node between brick nodes 6 and 24
3227  {
3228  unsigned j=15;
3229 
3230  // Need new node?
3231  Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
3232  old_node_pt=brick_edge_node_pt[edge];
3233  if (old_node_pt==0)
3234  {
3235  Node* new_node_pt=0;
3236  if (edge.is_boundary_edge())
3237  {
3238  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3239  }
3240  else
3241  {
3242  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3243  }
3244  brick_edge_node_pt[edge]=new_node_pt;
3245  Node_pt.push_back(new_node_pt);
3246  Vector<double> s(3);
3247  Vector<double> s_tet(3);
3248  Vector<double> x_tet(3);
3249  el_pt->local_coordinate_of_node(j,s);
3250  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3251  tet_el_pt->interpolated_x(s_tet,x_tet);
3252  new_node_pt->x(0)=x_tet[0];
3253  new_node_pt->x(1)=x_tet[1];
3254  new_node_pt->x(2)=x_tet[2];
3255  }
3256  // Node already exists
3257  else
3258  {
3259  el_pt->node_pt(j)=old_node_pt;
3260  }
3261  }
3262 
3263 
3264  // Brick edge node between brick nodes 8 and 26
3265  {
3266  unsigned j=17;
3267 
3268  // Need new node?
3269  Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
3270  old_node_pt=brick_edge_node_pt[edge];
3271  if (old_node_pt==0)
3272  {
3273  Node* new_node_pt=0;
3274  if (edge.is_boundary_edge())
3275  {
3276  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3277  }
3278  else
3279  {
3280  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3281  }
3282  brick_edge_node_pt[edge]=new_node_pt;
3283  Node_pt.push_back(new_node_pt);
3284  Vector<double> s(3);
3285  Vector<double> s_tet(3);
3286  Vector<double> x_tet(3);
3287  el_pt->local_coordinate_of_node(j,s);
3288  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3289  tet_el_pt->interpolated_x(s_tet,x_tet);
3290  new_node_pt->x(0)=x_tet[0];
3291  new_node_pt->x(1)=x_tet[1];
3292  new_node_pt->x(2)=x_tet[2];
3293  }
3294  // Node already exists
3295  else
3296  {
3297  el_pt->node_pt(j)=old_node_pt;
3298  }
3299  }
3300 
3301 
3302  // Mid brick-face node associated with face
3303  // spanned by mid-vertex nodes associated with tet edges 0 and 2
3304  {
3305  unsigned j=10;
3306 
3307  // Need new node?
3308  TFace face(tet_el_pt->node_pt( central_tet_vertex),
3309  tet_el_pt->node_pt(tet_edge_node[0][0]),
3310  tet_el_pt->node_pt(tet_edge_node[2][0]));
3311 
3312  old_node_pt=tet_face_node_pt[face];
3313  if (old_node_pt==0)
3314  {
3315  Node* new_node_pt=0;
3316  if (face.is_boundary_face())
3317  {
3318  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3319  }
3320  else
3321  {
3322  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3323  }
3324  tet_face_node_pt[face]=new_node_pt;
3325  Node_pt.push_back(new_node_pt);
3326  Vector<double> s(3);
3327  Vector<double> s_tet(3);
3328  Vector<double> x_tet(3);
3329  el_pt->local_coordinate_of_node(j,s);
3330  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3331  tet_el_pt->interpolated_x(s_tet,x_tet);
3332  new_node_pt->x(0)=x_tet[0];
3333  new_node_pt->x(1)=x_tet[1];
3334  new_node_pt->x(2)=x_tet[2];
3335  }
3336  // Node already exists
3337  else
3338  {
3339  el_pt->node_pt(j)=old_node_pt;
3340  }
3341  }
3342 
3343 
3344 
3345  // Mid brick-face node associated with face
3346  // spanned by mid-vertex nodes associated with tet edges 1 and 2
3347  {
3348  unsigned j=12;
3349 
3350  // Need new node?
3351  TFace face(tet_el_pt->node_pt(central_tet_vertex),
3352  tet_el_pt->node_pt(tet_edge_node[1][0]),
3353  tet_el_pt->node_pt(tet_edge_node[2][0]));
3354  old_node_pt=tet_face_node_pt[face];
3355  if (old_node_pt==0)
3356  {
3357  Node* new_node_pt=0;
3358  if (face.is_boundary_face())
3359  {
3360  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3361  }
3362  else
3363  {
3364  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3365  }
3366  tet_face_node_pt[face]=new_node_pt;
3367  Node_pt.push_back(new_node_pt);
3368  Vector<double> s(3);
3369  Vector<double> s_tet(3);
3370  Vector<double> x_tet(3);
3371  el_pt->local_coordinate_of_node(j,s);
3372  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3373  tet_el_pt->interpolated_x(s_tet,x_tet);
3374  new_node_pt->x(0)=x_tet[0];
3375  new_node_pt->x(1)=x_tet[1];
3376  new_node_pt->x(2)=x_tet[2];
3377  }
3378  // Node already exists
3379  else
3380  {
3381  el_pt->node_pt(j)=old_node_pt;
3382  }
3383  }
3384 
3385 
3386 
3387  // Mid brick-face node associated with face
3388  // spanned by mid-vertex nodes associated with tet edges 0 and 1
3389  {
3390  unsigned j=4;
3391 
3392  // Need new node?
3393  TFace face(tet_el_pt->node_pt(central_tet_vertex),
3394  tet_el_pt->node_pt(tet_edge_node[0][0]),
3395  tet_el_pt->node_pt(tet_edge_node[1][0]));
3396 
3397  old_node_pt=tet_face_node_pt[face];
3398  if (old_node_pt==0)
3399  {
3400  Node* new_node_pt=0;
3401  if (face.is_boundary_face())
3402  {
3403  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3404  }
3405  else
3406  {
3407  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3408  }
3409  tet_face_node_pt[face]=new_node_pt;
3410  Node_pt.push_back(new_node_pt);
3411  Vector<double> s(3);
3412  Vector<double> s_tet(3);
3413  Vector<double> x_tet(3);
3414  el_pt->local_coordinate_of_node(j,s);
3415  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3416  tet_el_pt->interpolated_x(s_tet,x_tet);
3417  new_node_pt->x(0)=x_tet[0];
3418  new_node_pt->x(1)=x_tet[1];
3419  new_node_pt->x(2)=x_tet[2];
3420  }
3421  // Node already exists
3422  else
3423  {
3424  el_pt->node_pt(j)=old_node_pt;
3425  }
3426  }
3427 
3428 
3429  // Top mid brick-face node -- only built by this element
3430  {
3431  unsigned j=22;
3432  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3433  Node_pt.push_back(new_node_pt);
3434  Vector<double> s(3);
3435  Vector<double> s_tet(3);
3436  Vector<double> x_tet(3);
3437  el_pt->local_coordinate_of_node(j,s);
3438  dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3439  top_mid_face_node2_pt=new_node_pt;
3440  tet_el_pt->interpolated_x(s_tet,x_tet);
3441  new_node_pt->x(0)=x_tet[0];
3442  new_node_pt->x(1)=x_tet[1];
3443  new_node_pt->x(2)=x_tet[2];
3444  }
3445 
3446 
3447  // Right mid brick-face node copy from first element
3448  {
3449  unsigned j=14;
3450 
3451  // Always copied
3452  el_pt->node_pt(j)=back_mid_face_node0_pt;
3453  }
3454 
3455 
3456  // Back mid brick-face node copy from previous element
3457  {
3458  unsigned j=16;
3459 
3460  // Always copied
3461  el_pt->node_pt(j)=right_mid_face_node1_pt;
3462  }
3463 
3464  }
3465 
3466 
3467  // Fourth brick element is centred at node 2 of tet:
3468  //--------------------------------------------------
3469  {
3470  // Assign coordinates of dummy element
3471  for (unsigned j=0;j<8;j++)
3472  {
3473  Node* nod_pt=dummy_q_el_pt[3]->node_pt(j);
3474  Vector<double> s_tet(3);
3475  Vector<double> x_tet(3);
3476  switch (j)
3477  {
3478  case 0:
3479  tet_el_pt->local_coordinate_of_node(2,s_tet);
3480  nod_pt->set_value(0,s_tet[0]);
3481  nod_pt->set_value(1,s_tet[1]);
3482  nod_pt->set_value(2,s_tet[2]);
3483  tet_el_pt->interpolated_x(s_tet,x_tet);
3484  nod_pt->x(0)=x_tet[0];
3485  nod_pt->x(1)=x_tet[1];
3486  nod_pt->x(2)=x_tet[2];
3487  break;
3488  case 1:
3489  tet_el_pt->local_coordinate_of_node(7,s_tet);
3490  nod_pt->set_value(0,s_tet[0]);
3491  nod_pt->set_value(1,s_tet[1]);
3492  nod_pt->set_value(2,s_tet[2]);
3493  tet_el_pt->interpolated_x(s_tet,x_tet);
3494  nod_pt->x(0)=x_tet[0];
3495  nod_pt->x(1)=x_tet[1];
3496  nod_pt->x(2)=x_tet[2];
3497  break;
3498  case 2:
3499  tet_el_pt->local_coordinate_of_node(5,s_tet);
3500  nod_pt->set_value(0,s_tet[0]);
3501  nod_pt->set_value(1,s_tet[1]);
3502  nod_pt->set_value(2,s_tet[2]);
3503  tet_el_pt->interpolated_x(s_tet,x_tet);
3504  nod_pt->x(0)=x_tet[0];
3505  nod_pt->x(1)=x_tet[1];
3506  nod_pt->x(2)=x_tet[2];
3507  break;
3508  case 3:
3509  // label 11 in initial sketch: Mid face node on face
3510  // spanned by tet nodes 0,1,2
3511  s_tet[0]=1.0/3.0;
3512  s_tet[1]=1.0/3.0;
3513  s_tet[2]=1.0/3.0;
3514  nod_pt->set_value(0,s_tet[0]);
3515  nod_pt->set_value(1,s_tet[1]);
3516  nod_pt->set_value(2,s_tet[2]);
3517  tet_el_pt->interpolated_x(s_tet,x_tet);
3518  nod_pt->x(0)=x_tet[0];
3519  nod_pt->x(1)=x_tet[1];
3520  nod_pt->x(2)=x_tet[2];
3521  break;
3522  case 4:
3523  tet_el_pt->local_coordinate_of_node(8,s_tet);
3524  nod_pt->set_value(0,s_tet[0]);
3525  nod_pt->set_value(1,s_tet[1]);
3526  nod_pt->set_value(2,s_tet[2]);
3527  tet_el_pt->interpolated_x(s_tet,x_tet);
3528  nod_pt->x(0)=x_tet[0];
3529  nod_pt->x(1)=x_tet[1];
3530  nod_pt->x(2)=x_tet[2];
3531  break;
3532  case 5:
3533  // label 10 in initial sketch: Mid face node on face
3534  // spanned by tet nodes 0,2,3
3535  s_tet[0]=0.0;
3536  s_tet[1]=1.0/3.0;
3537  s_tet[2]=1.0/3.0;
3538  nod_pt->set_value(0,s_tet[0]);
3539  nod_pt->set_value(1,s_tet[1]);
3540  nod_pt->set_value(2,s_tet[2]);
3541  tet_el_pt->interpolated_x(s_tet,x_tet);
3542  nod_pt->x(0)=x_tet[0];
3543  nod_pt->x(1)=x_tet[1];
3544  nod_pt->x(2)=x_tet[2];
3545  break;
3546  case 6:
3547  // label 12 in initial sketch: Mid face node on face
3548  // spanned by tet nodes 0,2,3
3549  s_tet[0]=1.0/3.0;
3550  s_tet[1]=0.0;
3551  s_tet[2]=1.0/3.0;
3552  nod_pt->set_value(0,s_tet[0]);
3553  nod_pt->set_value(1,s_tet[1]);
3554  nod_pt->set_value(2,s_tet[2]);
3555  tet_el_pt->interpolated_x(s_tet,x_tet);
3556  nod_pt->x(0)=x_tet[0];
3557  nod_pt->x(1)=x_tet[1];
3558  nod_pt->x(2)=x_tet[2];
3559  break;
3560  case 7:
3561  // label 14 in initial sketch: Centroid
3562  s_tet[0]=0.25;
3563  s_tet[1]=0.25;
3564  s_tet[2]=0.25;
3565  nod_pt->set_value(0,s_tet[0]);
3566  nod_pt->set_value(1,s_tet[1]);
3567  nod_pt->set_value(2,s_tet[2]);
3568  tet_el_pt->interpolated_x(s_tet,x_tet);
3569  nod_pt->x(0)=x_tet[0];
3570  nod_pt->x(1)=x_tet[1];
3571  nod_pt->x(2)=x_tet[2];
3572  break;
3573  }
3574  }
3575 
3576 
3577 
3578  // Create actual third brick element
3579  FiniteElement* el_pt=new ELEMENT;
3580  brick_el3_pt=el_pt;
3581  Element_pt.push_back(el_pt);
3582 
3583  TFace face0(tet_el_pt->node_pt(1),
3584  tet_el_pt->node_pt(2),
3585  tet_el_pt->node_pt(3));
3586 
3587  TFace face1(tet_el_pt->node_pt(0),
3588  tet_el_pt->node_pt(2),
3589  tet_el_pt->node_pt(3));
3590 
3591  TFace face2(tet_el_pt->node_pt(0),
3592  tet_el_pt->node_pt(1),
3593  tet_el_pt->node_pt(2));
3594 
3595  // Tet vertex nodes along edges emanating from node 0 in brick
3596  Vector<Vector<unsigned> > tet_edge_node(3);
3597  tet_edge_node[0].resize(2);
3598  tet_edge_node[0][0]=7;
3599  tet_edge_node[0][1]=1;
3600  tet_edge_node[1].resize(2);
3601  tet_edge_node[1][0]=5;
3602  tet_edge_node[1][1]=0;
3603  tet_edge_node[2].resize(2);
3604  tet_edge_node[2][0]=8;
3605  tet_edge_node[2][1]=3;
3606 
3607  // Node number of tet vertex that node 0 in brick is centred on
3608  unsigned central_tet_vertex=2;
3609 
3610  Node* tet_node_pt=0;
3611  Node* old_node_pt=0;
3612 
3613  // Corner node
3614  {
3615  unsigned j=0;
3616 
3617  // Need new node?
3618  tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
3619  old_node_pt=tet_node_node_pt[tet_node_pt];
3620  if (old_node_pt==0)
3621  {
3622  Node* new_node_pt=0;
3623  if (tet_node_pt->is_on_boundary())
3624  {
3625  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3626  }
3627  else
3628  {
3629  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3630  }
3631  tet_node_node_pt[tet_node_pt]=new_node_pt;
3632  Node_pt.push_back(new_node_pt);
3633  Vector<double> s(3);
3634  Vector<double> s_tet(3);
3635  Vector<double> x_tet(3);
3636  el_pt->local_coordinate_of_node(j,s);
3637  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3638  tet_el_pt->interpolated_x(s_tet,x_tet);
3639  new_node_pt->x(0)=x_tet[0];
3640  new_node_pt->x(1)=x_tet[1];
3641  new_node_pt->x(2)=x_tet[2];
3642  }
3643  // Node already exists
3644  else
3645  {
3646  el_pt->node_pt(j)=old_node_pt;
3647  }
3648  }
3649 
3650 
3651  // Brick vertex node coindides with mid-edge node on tet edge 0
3652  {
3653  unsigned j=2;
3654 
3655  // Need new node?
3656  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
3657  old_node_pt=tet_node_node_pt[tet_node_pt];
3658  if (old_node_pt==0)
3659  {
3660  Node* new_node_pt=0;
3661  if (tet_node_pt->is_on_boundary())
3662  {
3663  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3664  }
3665  else
3666  {
3667  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3668  }
3669  tet_node_node_pt[tet_node_pt]=new_node_pt;
3670  Node_pt.push_back(new_node_pt);
3671  Vector<double> s(3);
3672  Vector<double> s_tet(3);
3673  Vector<double> x_tet(3);
3674  el_pt->local_coordinate_of_node(j,s);
3675  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3676  tet_el_pt->interpolated_x(s_tet,x_tet);
3677  new_node_pt->x(0)=x_tet[0];
3678  new_node_pt->x(1)=x_tet[1];
3679  new_node_pt->x(2)=x_tet[2];
3680  }
3681  // Node already exists
3682  else
3683  {
3684  el_pt->node_pt(j)=old_node_pt;
3685  }
3686  }
3687 
3688 
3689  // Brick vertex node coindides with mid vertex node of tet edge 1
3690  {
3691  unsigned j=6;
3692 
3693  // Need new node?
3694  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
3695  old_node_pt=tet_node_node_pt[tet_node_pt];
3696  if (old_node_pt==0)
3697  {
3698  Node* new_node_pt=0;
3699  if (tet_node_pt->is_on_boundary())
3700  {
3701  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3702  }
3703  else
3704  {
3705  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3706  }
3707  tet_node_node_pt[tet_node_pt]=new_node_pt;
3708  Node_pt.push_back(new_node_pt);
3709  Vector<double> s(3);
3710  Vector<double> s_tet(3);
3711  Vector<double> x_tet(3);
3712  el_pt->local_coordinate_of_node(j,s);
3713  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3714  tet_el_pt->interpolated_x(s_tet,x_tet);
3715  new_node_pt->x(0)=x_tet[0];
3716  new_node_pt->x(1)=x_tet[1];
3717  new_node_pt->x(2)=x_tet[2];
3718  }
3719  // Node already exists
3720  else
3721  {
3722  el_pt->node_pt(j)=old_node_pt;
3723  }
3724  }
3725 
3726 
3727  // Brick vertex node coindides with mid-vertex node of tet edge 2
3728  {
3729  unsigned j=18;
3730 
3731  // Need new node?
3732  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
3733  old_node_pt=tet_node_node_pt[tet_node_pt];
3734  if (old_node_pt==0)
3735  {
3736  Node* new_node_pt=0;
3737  if (tet_node_pt->is_on_boundary())
3738  {
3739  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3740  }
3741  else
3742  {
3743  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3744  }
3745  tet_node_node_pt[tet_node_pt]=new_node_pt;
3746  Node_pt.push_back(new_node_pt);
3747  Vector<double> s(3);
3748  Vector<double> s_tet(3);
3749  Vector<double> x_tet(3);
3750  el_pt->local_coordinate_of_node(j,s);
3751  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3752  tet_el_pt->interpolated_x(s_tet,x_tet);
3753  new_node_pt->x(0)=x_tet[0];
3754  new_node_pt->x(1)=x_tet[1];
3755  new_node_pt->x(2)=x_tet[2];
3756  }
3757  // Node already exists
3758  else
3759  {
3760  el_pt->node_pt(j)=old_node_pt;
3761  }
3762  }
3763 
3764 
3765 
3766  // Brick vertex node in the middle of tet face0
3767  {
3768  unsigned j=20;
3769 
3770  // Need new node?
3771  old_node_pt=tet_face_node_pt[face0];
3772  if (old_node_pt==0)
3773  {
3774  Node* new_node_pt=0;
3775  if (face0.is_boundary_face())
3776  {
3777  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3778  }
3779  else
3780  {
3781  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3782  }
3783  tet_face_node_pt[face0]=new_node_pt;
3784  Node_pt.push_back(new_node_pt);
3785  Vector<double> s(3);
3786  Vector<double> s_tet(3);
3787  Vector<double> x_tet(3);
3788  el_pt->local_coordinate_of_node(j,s);
3789  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3790  tet_el_pt->interpolated_x(s_tet,x_tet);
3791  new_node_pt->x(0)=x_tet[0];
3792  new_node_pt->x(1)=x_tet[1];
3793  new_node_pt->x(2)=x_tet[2];
3794  }
3795  // Node already exists
3796  else
3797  {
3798  el_pt->node_pt(j)=old_node_pt;
3799  }
3800  }
3801 
3802  // Brick vertex node in the middle of tet face1
3803  {
3804  unsigned j=24;
3805 
3806  // Need new node?
3807  old_node_pt=tet_face_node_pt[face1];
3808  if (old_node_pt==0)
3809  {
3810  Node* new_node_pt=0;
3811  if (face1.is_boundary_face())
3812  {
3813  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3814  }
3815  else
3816  {
3817  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3818  }
3819  tet_face_node_pt[face1]=new_node_pt;
3820  Node_pt.push_back(new_node_pt);
3821  Vector<double> s(3);
3822  Vector<double> s_tet(3);
3823  Vector<double> x_tet(3);
3824  el_pt->local_coordinate_of_node(j,s);
3825  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3826  tet_el_pt->interpolated_x(s_tet,x_tet);
3827  new_node_pt->x(0)=x_tet[0];
3828  new_node_pt->x(1)=x_tet[1];
3829  new_node_pt->x(2)=x_tet[2];
3830  }
3831  // Node already exists
3832  else
3833  {
3834  el_pt->node_pt(j)=old_node_pt;
3835  }
3836  }
3837 
3838  // Brick vertex node in the middle of tet face2
3839  {
3840  unsigned j=8;
3841 
3842  // Need new node?
3843  old_node_pt=tet_face_node_pt[face2];
3844  if (old_node_pt==0)
3845  {
3846  Node* new_node_pt=0;
3847  if (face2.is_boundary_face())
3848  {
3849  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3850  }
3851  else
3852  {
3853  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3854  }
3855  tet_face_node_pt[face2]=new_node_pt;
3856  Node_pt.push_back(new_node_pt);
3857  Vector<double> s(3);
3858  Vector<double> s_tet(3);
3859  Vector<double> x_tet(3);
3860  el_pt->local_coordinate_of_node(j,s);
3861  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3862  tet_el_pt->interpolated_x(s_tet,x_tet);
3863  new_node_pt->x(0)=x_tet[0];
3864  new_node_pt->x(1)=x_tet[1];
3865  new_node_pt->x(2)=x_tet[2];
3866  }
3867  // Node already exists
3868  else
3869  {
3870  el_pt->node_pt(j)=old_node_pt;
3871  }
3872  }
3873 
3874 
3875  // Brick vertex node in centroid of tet. Only built for first element.
3876  // Enumerated "13" in initial sketch.
3877  {
3878  unsigned j=26;
3879 
3880  // Always copied
3881  el_pt->node_pt(j)=centroid_node_pt;
3882  }
3883 
3884 
3885  // Internal brick node -- always built
3886  {
3887  unsigned j=13;
3888 
3889  // Always new
3890  {
3891  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3892  Node_pt.push_back(new_node_pt);
3893  Vector<double> s(3);
3894  Vector<double> s_tet(3);
3895  Vector<double> x_tet(3);
3896  el_pt->local_coordinate_of_node(j,s);
3897  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3898  tet_el_pt->interpolated_x(s_tet,x_tet);
3899  new_node_pt->x(0)=x_tet[0];
3900  new_node_pt->x(1)=x_tet[1];
3901  new_node_pt->x(2)=x_tet[2];
3902  }
3903  }
3904 
3905  // Brick edge node between brick nodes 0 and 2
3906  {
3907  unsigned j=1;
3908 
3909  // Need new node?
3910  Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
3911  old_node_pt=brick_edge_node_pt[edge];
3912  if (old_node_pt==0)
3913  {
3914  Node* new_node_pt=0;
3915  if (edge.is_boundary_edge())
3916  {
3917  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3918  }
3919  else
3920  {
3921  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3922  }
3923  brick_edge_node_pt[edge]=new_node_pt;
3924  Node_pt.push_back(new_node_pt);
3925  Vector<double> s(3);
3926  Vector<double> s_tet(3);
3927  Vector<double> x_tet(3);
3928  el_pt->local_coordinate_of_node(j,s);
3929  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3930  tet_el_pt->interpolated_x(s_tet,x_tet);
3931  new_node_pt->x(0)=x_tet[0];
3932  new_node_pt->x(1)=x_tet[1];
3933  new_node_pt->x(2)=x_tet[2];
3934  }
3935  // Node already exists
3936  else
3937  {
3938  el_pt->node_pt(j)=old_node_pt;
3939  }
3940  }
3941 
3942 
3943  // Brick edge node between brick nodes 0 and 6
3944  {
3945  unsigned j=3;
3946 
3947  // Need new node?
3948  Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
3949  old_node_pt=brick_edge_node_pt[edge];
3950  if (old_node_pt==0)
3951  {
3952  Node* new_node_pt=0;
3953  if (edge.is_boundary_edge())
3954  {
3955  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3956  }
3957  else
3958  {
3959  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3960  }
3961  brick_edge_node_pt[edge]=new_node_pt;
3962  Node_pt.push_back(new_node_pt);
3963  Vector<double> s(3);
3964  Vector<double> s_tet(3);
3965  Vector<double> x_tet(3);
3966  el_pt->local_coordinate_of_node(j,s);
3967  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3968  tet_el_pt->interpolated_x(s_tet,x_tet);
3969  new_node_pt->x(0)=x_tet[0];
3970  new_node_pt->x(1)=x_tet[1];
3971  new_node_pt->x(2)=x_tet[2];
3972  }
3973  // Node already exists
3974  else
3975  {
3976  el_pt->node_pt(j)=old_node_pt;
3977  }
3978  }
3979 
3980  // Brick edge node between brick nodes 2 and 8
3981  {
3982  unsigned j=5;
3983 
3984  // Need new node?
3985  Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
3986  old_node_pt=brick_edge_node_pt[edge];
3987  if (old_node_pt==0)
3988  {
3989  Node* new_node_pt=0;
3990  if (edge.is_boundary_edge())
3991  {
3992  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3993  }
3994  else
3995  {
3996  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3997  }
3998  brick_edge_node_pt[edge]=new_node_pt;
3999  Node_pt.push_back(new_node_pt);
4000  Vector<double> s(3);
4001  Vector<double> s_tet(3);
4002  Vector<double> x_tet(3);
4003  el_pt->local_coordinate_of_node(j,s);
4004  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4005  tet_el_pt->interpolated_x(s_tet,x_tet);
4006  new_node_pt->x(0)=x_tet[0];
4007  new_node_pt->x(1)=x_tet[1];
4008  new_node_pt->x(2)=x_tet[2];
4009  }
4010  // Node already exists
4011  else
4012  {
4013  el_pt->node_pt(j)=old_node_pt;
4014  }
4015  }
4016 
4017  // Brick edge node between brick nodes 6 and 8
4018  {
4019  unsigned j=7;
4020 
4021  // Need new node?
4022  Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
4023  old_node_pt=brick_edge_node_pt[edge];
4024  if (old_node_pt==0)
4025  {
4026  Node* new_node_pt=0;
4027  if (edge.is_boundary_edge())
4028  {
4029  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4030  }
4031  else
4032  {
4033  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4034  }
4035  brick_edge_node_pt[edge]=new_node_pt;
4036  Node_pt.push_back(new_node_pt);
4037  Vector<double> s(3);
4038  Vector<double> s_tet(3);
4039  Vector<double> x_tet(3);
4040  el_pt->local_coordinate_of_node(j,s);
4041  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4042  tet_el_pt->interpolated_x(s_tet,x_tet);
4043  new_node_pt->x(0)=x_tet[0];
4044  new_node_pt->x(1)=x_tet[1];
4045  new_node_pt->x(2)=x_tet[2];
4046  }
4047  // Node already exists
4048  else
4049  {
4050  el_pt->node_pt(j)=old_node_pt;
4051  }
4052  }
4053 
4054  // Brick edge node between brick nodes 18 and 20
4055  {
4056  unsigned j=19;
4057 
4058  // Need new node?
4059  Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
4060  old_node_pt=brick_edge_node_pt[edge];
4061  if (old_node_pt==0)
4062  {
4063  Node* new_node_pt=0;
4064  if (edge.is_boundary_edge())
4065  {
4066  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4067  }
4068  else
4069  {
4070  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4071  }
4072  brick_edge_node_pt[edge]=new_node_pt;
4073  Node_pt.push_back(new_node_pt);
4074  Vector<double> s(3);
4075  Vector<double> s_tet(3);
4076  Vector<double> x_tet(3);
4077  el_pt->local_coordinate_of_node(j,s);
4078  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4079  tet_el_pt->interpolated_x(s_tet,x_tet);
4080  new_node_pt->x(0)=x_tet[0];
4081  new_node_pt->x(1)=x_tet[1];
4082  new_node_pt->x(2)=x_tet[2];
4083  }
4084  // Node already exists
4085  else
4086  {
4087  el_pt->node_pt(j)=old_node_pt;
4088  }
4089  }
4090 
4091 
4092  // Brick edge node between brick nodes 18 and 24
4093  {
4094  unsigned j=21;
4095 
4096  // Need new node?
4097  Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
4098  old_node_pt=brick_edge_node_pt[edge];
4099  if (old_node_pt==0)
4100  {
4101  Node* new_node_pt=0;
4102  if (edge.is_boundary_edge())
4103  {
4104  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4105  }
4106  else
4107  {
4108  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4109  }
4110  brick_edge_node_pt[edge]=new_node_pt;
4111  Node_pt.push_back(new_node_pt);
4112  Vector<double> s(3);
4113  Vector<double> s_tet(3);
4114  Vector<double> x_tet(3);
4115  el_pt->local_coordinate_of_node(j,s);
4116  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4117  tet_el_pt->interpolated_x(s_tet,x_tet);
4118  new_node_pt->x(0)=x_tet[0];
4119  new_node_pt->x(1)=x_tet[1];
4120  new_node_pt->x(2)=x_tet[2];
4121  }
4122  // Node already exists
4123  else
4124  {
4125  el_pt->node_pt(j)=old_node_pt;
4126  }
4127  }
4128 
4129  // Brick edge node between brick nodes 20 and 26
4130  {
4131  unsigned j=23;
4132 
4133  // Need new node?
4134  Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
4135  old_node_pt=brick_edge_node_pt[edge];
4136  if (old_node_pt==0)
4137  {
4138  Node* new_node_pt=0;
4139  if (edge.is_boundary_edge())
4140  {
4141  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4142  }
4143  else
4144  {
4145  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4146  }
4147  brick_edge_node_pt[edge]=new_node_pt;
4148  Node_pt.push_back(new_node_pt);
4149  Vector<double> s(3);
4150  Vector<double> s_tet(3);
4151  Vector<double> x_tet(3);
4152  el_pt->local_coordinate_of_node(j,s);
4153  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4154  tet_el_pt->interpolated_x(s_tet,x_tet);
4155  new_node_pt->x(0)=x_tet[0];
4156  new_node_pt->x(1)=x_tet[1];
4157  new_node_pt->x(2)=x_tet[2];
4158  }
4159  // Node already exists
4160  else
4161  {
4162  el_pt->node_pt(j)=old_node_pt;
4163  }
4164  }
4165 
4166 
4167  // Brick edge node between brick nodes 24 and 26
4168  {
4169  unsigned j=25;
4170 
4171  // Need new node?
4172  Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
4173  old_node_pt=brick_edge_node_pt[edge];
4174  if (old_node_pt==0)
4175  {
4176  Node* new_node_pt=0;
4177  if (edge.is_boundary_edge())
4178  {
4179  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4180  }
4181  else
4182  {
4183  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4184  }
4185  brick_edge_node_pt[edge]=new_node_pt;
4186  Node_pt.push_back(new_node_pt);
4187  Vector<double> s(3);
4188  Vector<double> s_tet(3);
4189  Vector<double> x_tet(3);
4190  el_pt->local_coordinate_of_node(j,s);
4191  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4192  tet_el_pt->interpolated_x(s_tet,x_tet);
4193  new_node_pt->x(0)=x_tet[0];
4194  new_node_pt->x(1)=x_tet[1];
4195  new_node_pt->x(2)=x_tet[2];
4196  }
4197  // Node already exists
4198  else
4199  {
4200  el_pt->node_pt(j)=old_node_pt;
4201  }
4202  }
4203 
4204  // Brick edge node between brick nodes 0 and 18
4205  {
4206  unsigned j=9;
4207 
4208  // Need new node?
4209  Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
4210  old_node_pt=brick_edge_node_pt[edge];
4211  if (old_node_pt==0)
4212  {
4213  Node* new_node_pt=0;
4214  if (edge.is_boundary_edge())
4215  {
4216  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4217  }
4218  else
4219  {
4220  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4221  }
4222  brick_edge_node_pt[edge]=new_node_pt;
4223  Node_pt.push_back(new_node_pt);
4224  Vector<double> s(3);
4225  Vector<double> s_tet(3);
4226  Vector<double> x_tet(3);
4227  el_pt->local_coordinate_of_node(j,s);
4228  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4229  tet_el_pt->interpolated_x(s_tet,x_tet);
4230  new_node_pt->x(0)=x_tet[0];
4231  new_node_pt->x(1)=x_tet[1];
4232  new_node_pt->x(2)=x_tet[2];
4233  }
4234  // Node already exists
4235  else
4236  {
4237  el_pt->node_pt(j)=old_node_pt;
4238  }
4239  }
4240 
4241 
4242  // Brick edge node between brick nodes 2 and 20
4243  {
4244  unsigned j=11;
4245 
4246  // Need new node?
4247  Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
4248  old_node_pt=brick_edge_node_pt[edge];
4249  if (old_node_pt==0)
4250  {
4251  Node* new_node_pt=0;
4252  if (edge.is_boundary_edge())
4253  {
4254  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4255  }
4256  else
4257  {
4258  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4259  }
4260  brick_edge_node_pt[edge]=new_node_pt;
4261  Node_pt.push_back(new_node_pt);
4262  Vector<double> s(3);
4263  Vector<double> s_tet(3);
4264  Vector<double> x_tet(3);
4265  el_pt->local_coordinate_of_node(j,s);
4266  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4267  tet_el_pt->interpolated_x(s_tet,x_tet);
4268  new_node_pt->x(0)=x_tet[0];
4269  new_node_pt->x(1)=x_tet[1];
4270  new_node_pt->x(2)=x_tet[2];
4271  }
4272  // Node already exists
4273  else
4274  {
4275  el_pt->node_pt(j)=old_node_pt;
4276  }
4277  }
4278 
4279 
4280  // Brick edge node between brick nodes 6 and 24
4281  {
4282  unsigned j=15;
4283 
4284  // Need new node?
4285  Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
4286  old_node_pt=brick_edge_node_pt[edge];
4287  if (old_node_pt==0)
4288  {
4289  Node* new_node_pt=0;
4290  if (edge.is_boundary_edge())
4291  {
4292  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4293  }
4294  else
4295  {
4296  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4297  }
4298  brick_edge_node_pt[edge]=new_node_pt;
4299  Node_pt.push_back(new_node_pt);
4300  Vector<double> s(3);
4301  Vector<double> s_tet(3);
4302  Vector<double> x_tet(3);
4303  el_pt->local_coordinate_of_node(j,s);
4304  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4305  tet_el_pt->interpolated_x(s_tet,x_tet);
4306  new_node_pt->x(0)=x_tet[0];
4307  new_node_pt->x(1)=x_tet[1];
4308  new_node_pt->x(2)=x_tet[2];
4309  }
4310  // Node already exists
4311  else
4312  {
4313  el_pt->node_pt(j)=old_node_pt;
4314  }
4315  }
4316 
4317 
4318  // Brick edge node between brick nodes 8 and 26
4319  {
4320  unsigned j=17;
4321 
4322  // Need new node?
4323  Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
4324  old_node_pt=brick_edge_node_pt[edge];
4325  if (old_node_pt==0)
4326  {
4327  Node* new_node_pt=0;
4328  if (edge.is_boundary_edge())
4329  {
4330  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4331  }
4332  else
4333  {
4334  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4335  }
4336  brick_edge_node_pt[edge]=new_node_pt;
4337  Node_pt.push_back(new_node_pt);
4338  Vector<double> s(3);
4339  Vector<double> s_tet(3);
4340  Vector<double> x_tet(3);
4341  el_pt->local_coordinate_of_node(j,s);
4342  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4343  tet_el_pt->interpolated_x(s_tet,x_tet);
4344  new_node_pt->x(0)=x_tet[0];
4345  new_node_pt->x(1)=x_tet[1];
4346  new_node_pt->x(2)=x_tet[2];
4347  }
4348  // Node already exists
4349  else
4350  {
4351  el_pt->node_pt(j)=old_node_pt;
4352  }
4353  }
4354 
4355 
4356  // Mid brick-face node associated with face
4357  // spanned by mid-vertex nodes associated with tet edges 0 and 2
4358  {
4359  unsigned j=10;
4360 
4361  // Need new node?
4362  TFace face(tet_el_pt->node_pt( central_tet_vertex),
4363  tet_el_pt->node_pt(tet_edge_node[0][0]),
4364  tet_el_pt->node_pt(tet_edge_node[2][0]));
4365 
4366  old_node_pt=tet_face_node_pt[face];
4367  if (old_node_pt==0)
4368  {
4369  Node* new_node_pt=0;
4370  if (face.is_boundary_face())
4371  {
4372  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4373  }
4374  else
4375  {
4376  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4377  }
4378  tet_face_node_pt[face]=new_node_pt;
4379  Node_pt.push_back(new_node_pt);
4380  Vector<double> s(3);
4381  Vector<double> s_tet(3);
4382  Vector<double> x_tet(3);
4383  el_pt->local_coordinate_of_node(j,s);
4384  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4385  tet_el_pt->interpolated_x(s_tet,x_tet);
4386  new_node_pt->x(0)=x_tet[0];
4387  new_node_pt->x(1)=x_tet[1];
4388  new_node_pt->x(2)=x_tet[2];
4389  }
4390  // Node already exists
4391  else
4392  {
4393  el_pt->node_pt(j)=old_node_pt;
4394  }
4395  }
4396 
4397 
4398 
4399  // Mid brick-face node associated with face
4400  // spanned by mid-vertex nodes associated with tet edges 1 and 2
4401  {
4402  unsigned j=12;
4403 
4404  // Need new node?
4405  TFace face(tet_el_pt->node_pt(central_tet_vertex),
4406  tet_el_pt->node_pt(tet_edge_node[1][0]),
4407  tet_el_pt->node_pt(tet_edge_node[2][0]));
4408 
4409  old_node_pt=tet_face_node_pt[face];
4410  if (old_node_pt==0)
4411  {
4412  Node* new_node_pt=0;
4413  if (face.is_boundary_face())
4414  {
4415  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4416  }
4417  else
4418  {
4419  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4420  }
4421  tet_face_node_pt[face]=new_node_pt;
4422  Node_pt.push_back(new_node_pt);
4423  Vector<double> s(3);
4424  Vector<double> s_tet(3);
4425  Vector<double> x_tet(3);
4426  el_pt->local_coordinate_of_node(j,s);
4427  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4428  tet_el_pt->interpolated_x(s_tet,x_tet);
4429  new_node_pt->x(0)=x_tet[0];
4430  new_node_pt->x(1)=x_tet[1];
4431  new_node_pt->x(2)=x_tet[2];
4432  }
4433  // Node already exists
4434  else
4435  {
4436  el_pt->node_pt(j)=old_node_pt;
4437  }
4438  }
4439 
4440 
4441 
4442  // Mid brick-face node associated with face
4443  // spanned by mid-vertex nodes associated with tet edges 0 and 1
4444  {
4445  unsigned j=4;
4446 
4447  // Need new node?
4448  TFace face(tet_el_pt->node_pt(central_tet_vertex),
4449  tet_el_pt->node_pt(tet_edge_node[0][0]),
4450  tet_el_pt->node_pt(tet_edge_node[1][0]));
4451 
4452  old_node_pt=tet_face_node_pt[face];
4453  if (old_node_pt==0)
4454  {
4455  Node* new_node_pt=0;
4456  if (face.is_boundary_face())
4457  {
4458  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4459  }
4460  else
4461  {
4462  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4463  }
4464  tet_face_node_pt[face]=new_node_pt;
4465  Node_pt.push_back(new_node_pt);
4466  Vector<double> s(3);
4467  Vector<double> s_tet(3);
4468  Vector<double> x_tet(3);
4469  el_pt->local_coordinate_of_node(j,s);
4470  dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4471  tet_el_pt->interpolated_x(s_tet,x_tet);
4472  new_node_pt->x(0)=x_tet[0];
4473  new_node_pt->x(1)=x_tet[1];
4474  new_node_pt->x(2)=x_tet[2];
4475  }
4476  // Node already exists
4477  else
4478  {
4479  el_pt->node_pt(j)=old_node_pt;
4480  }
4481  }
4482 
4483 
4484  // Top mid brick-face node copy from top of second element
4485  {
4486  unsigned j=22;
4487 
4488  // Always copied
4489  el_pt->node_pt(j)=top_mid_face_node2_pt;
4490  }
4491 
4492 
4493  // Right mid brick-face node copy from top of first element
4494  {
4495  unsigned j=14;
4496 
4497  // Always copied
4498  el_pt->node_pt(j)=top_mid_face_node1_pt;
4499  }
4500 
4501 
4502  // Back mid brick-face node copy from top of zeroth element
4503  {
4504  unsigned j=16;
4505 
4506  // Always copied
4507  el_pt->node_pt(j)=top_mid_face_node0_pt;
4508  }
4509 
4510  }
4511 
4512 
4513 
4514  // Check if the four faces of the tet are on a boundary.
4515  // If they are, add the nodes in the brick mesh to those
4516  // boundaries.
4517 
4518  // Loop over faces of tet (in its own face index enumeration)
4519  for (int face_index=0;face_index<4;face_index++)
4520  {
4521 
4522  // Identify face and coordinates in it
4523  TFace* face_pt=0;
4524  switch (face_index)
4525  {
4526 
4527  case 0:
4528  // Face 0: s[0]=0; opposite node 0
4529  face_pt=new TFace(tet_el_pt->node_pt(1),
4530  tet_el_pt->node_pt(2),
4531  tet_el_pt->node_pt(3));
4532  break;
4533 
4534  case 1:
4535  // Face 1: s[1]=0; opposite node 1
4536  face_pt=new TFace(tet_el_pt->node_pt(0),
4537  tet_el_pt->node_pt(2),
4538  tet_el_pt->node_pt(3));
4539 
4540  break;
4541 
4542 
4543  case 2:
4544  // Face 2: s[2]=0; opposite node 2
4545  face_pt=new TFace(tet_el_pt->node_pt(0),
4546  tet_el_pt->node_pt(1),
4547  tet_el_pt->node_pt(3));
4548 
4549  break;
4550 
4551  case 3:
4552  // Face 3: s[0]+s[1]+s[2]=1; opposite node 3
4553  face_pt=new TFace(tet_el_pt->node_pt(0),
4554  tet_el_pt->node_pt(1),
4555  tet_el_pt->node_pt(2));
4556  break;
4557 
4558  }
4559 
4560 
4561  if (face_pt->is_boundary_face())
4562  {
4563  std::set<unsigned> bnd;
4564  std::set<unsigned>* bnd_pt=&bnd;
4565  face_pt->get_boundaries_pt(bnd_pt);
4566 
4567 #ifdef PARANOID
4568  if ((*bnd_pt).size()>1)
4569  {
4570  std::ostringstream error_stream;
4571  error_stream << "TFace should only be on one boundary.\n";
4572  throw OomphLibError(
4573  error_stream.str(),
4574  OOMPH_CURRENT_FUNCTION,
4575  OOMPH_EXCEPTION_LOCATION);
4576  }
4577 #endif
4578 
4579  if ((*bnd_pt).size()==1)
4580  {
4581  // Create new face element
4582  FaceElement* face_el_pt=0;
4583  if (tet_mesh_is_solid_mesh)
4584  {
4585 #ifdef PARANOID
4586  if (dynamic_cast<SolidTElement<3,3>*>(tet_el_pt)==0)
4587  {
4588  std::ostringstream error_stream;
4589  error_stream
4590  << "Tet-element cannot be cast to SolidTElement<3,3>.\n"
4591  << "BrickFromTetMesh can only be built from\n"
4592  << "mesh containing quadratic tets.\n"
4593  << std::endl;
4594  throw OomphLibError(error_stream.str(),
4595  OOMPH_CURRENT_FUNCTION,
4596  OOMPH_EXCEPTION_LOCATION);
4597  }
4598 #endif
4599  // Build the face element
4600  face_el_pt=
4601  new DummyFaceElement<SolidTElement<3,3> >(tet_el_pt,face_index);
4602  }
4603  else
4604  {
4605 #ifdef PARANOID
4606  if (dynamic_cast<TElement<3,3>*>(tet_el_pt)==0)
4607  {
4608  std::ostringstream error_stream;
4609  error_stream
4610  << "Tet-element cannot be cast to TElement<3,3>.\n"
4611  << "BrickFromTetMesh can only be built from\n"
4612  << "mesh containing quadratic tets.\n"
4613  << std::endl;
4614  throw OomphLibError(error_stream.str(),
4615  OOMPH_CURRENT_FUNCTION,
4616  OOMPH_EXCEPTION_LOCATION);
4617  }
4618 #endif
4619 
4620  // Build the face element
4621  face_el_pt=
4622  new DummyFaceElement<TElement<3,3> >(tet_el_pt,face_index);
4623  }
4624 
4625 
4626  // Specify boundary id in bulk mesh (needed to extract
4627  // boundary coordinate)
4628  unsigned b=(*(*bnd_pt).begin());
4629  Boundary_coordinate_exists[b]=true;
4630  face_el_pt->set_boundary_number_in_bulk_mesh(b);
4631 
4632  // Now set up the brick nodes on this face, enumerated as
4633  // in s_face
4634  Vector<Node*> brick_face_node_pt(19);
4635 
4636  switch (face_index)
4637  {
4638  case 0:
4639  brick_face_node_pt[0]=brick_el1_pt->node_pt(0);
4640  brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
4641  brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
4642 
4643  brick_face_node_pt[3]=brick_el1_pt->node_pt(18);
4644  brick_face_node_pt[4]=brick_el2_pt->node_pt(18);
4645  brick_face_node_pt[5]=brick_el1_pt->node_pt(2);
4646 
4647  brick_face_node_pt[6]=brick_el1_pt->node_pt(9);
4648  brick_face_node_pt[7]=brick_el3_pt->node_pt(1);
4649  brick_face_node_pt[8]=brick_el3_pt->node_pt(9);
4650 
4651  brick_face_node_pt[9]=brick_el2_pt->node_pt(9);
4652  brick_face_node_pt[10]=brick_el2_pt->node_pt(3);
4653  brick_face_node_pt[11]=brick_el1_pt->node_pt(1);
4654 
4655  brick_face_node_pt[12]=brick_el1_pt->node_pt(20);
4656 
4657  brick_face_node_pt[13]=brick_el2_pt->node_pt(12);
4658  brick_face_node_pt[14]=brick_el1_pt->node_pt(19);
4659 
4660  brick_face_node_pt[15]=brick_el1_pt->node_pt(10);
4661  brick_face_node_pt[16]=brick_el2_pt->node_pt(21);
4662 
4663  brick_face_node_pt[17]=brick_el3_pt->node_pt(10);
4664  brick_face_node_pt[18]=brick_el1_pt->node_pt(11);
4665  break;
4666 
4667  case 1:
4668  brick_face_node_pt[0]=brick_el0_pt->node_pt(0);
4669  brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
4670  brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
4671 
4672  brick_face_node_pt[3]=brick_el0_pt->node_pt(18);
4673  brick_face_node_pt[4]=brick_el2_pt->node_pt(18);
4674  brick_face_node_pt[5]=brick_el0_pt->node_pt(6);
4675 
4676  brick_face_node_pt[6]=brick_el0_pt->node_pt(9);
4677  brick_face_node_pt[7]=brick_el3_pt->node_pt(3);
4678  brick_face_node_pt[8]=brick_el3_pt->node_pt(9);
4679 
4680  brick_face_node_pt[9]=brick_el2_pt->node_pt(9);
4681  brick_face_node_pt[10]=brick_el2_pt->node_pt(1);
4682  brick_face_node_pt[11]=brick_el0_pt->node_pt(3);
4683 
4684  brick_face_node_pt[12]=brick_el0_pt->node_pt(24);
4685 
4686  brick_face_node_pt[13]=brick_el2_pt->node_pt(10);
4687  brick_face_node_pt[14]=brick_el0_pt->node_pt(21);
4688 
4689  brick_face_node_pt[15]=brick_el0_pt->node_pt(12);
4690  brick_face_node_pt[16]=brick_el3_pt->node_pt(21);
4691 
4692  brick_face_node_pt[17]=brick_el3_pt->node_pt(12);
4693  brick_face_node_pt[18]=brick_el0_pt->node_pt(15);
4694  break;
4695 
4696  case 2:
4697  brick_face_node_pt[0]=brick_el0_pt->node_pt(0);
4698  brick_face_node_pt[1]=brick_el1_pt->node_pt(0);
4699  brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
4700 
4701  brick_face_node_pt[3]=brick_el0_pt->node_pt(2);
4702  brick_face_node_pt[4]=brick_el1_pt->node_pt(2);
4703  brick_face_node_pt[5]=brick_el0_pt->node_pt(6);
4704 
4705  brick_face_node_pt[6]=brick_el0_pt->node_pt(1);
4706  brick_face_node_pt[7]=brick_el1_pt->node_pt(3);
4707  brick_face_node_pt[8]=brick_el1_pt->node_pt(1);
4708 
4709  brick_face_node_pt[9]=brick_el2_pt->node_pt(3);
4710  brick_face_node_pt[10]=brick_el2_pt->node_pt(1);
4711  brick_face_node_pt[11]=brick_el0_pt->node_pt(3);
4712 
4713  brick_face_node_pt[12]=brick_el0_pt->node_pt(8);
4714 
4715  brick_face_node_pt[13]=brick_el2_pt->node_pt(4);
4716  brick_face_node_pt[14]=brick_el0_pt->node_pt(5);
4717 
4718  brick_face_node_pt[15]=brick_el0_pt->node_pt(4);
4719  brick_face_node_pt[16]=brick_el1_pt->node_pt(5);
4720 
4721  brick_face_node_pt[17]=brick_el1_pt->node_pt(4);
4722  brick_face_node_pt[18]=brick_el0_pt->node_pt(7);
4723  break;
4724 
4725  case 3:
4726  brick_face_node_pt[0]=brick_el1_pt->node_pt(0);
4727  brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
4728  brick_face_node_pt[2]=brick_el0_pt->node_pt(0);
4729 
4730  brick_face_node_pt[3]=brick_el1_pt->node_pt(18);
4731  brick_face_node_pt[4]=brick_el3_pt->node_pt(6);
4732  brick_face_node_pt[5]=brick_el1_pt->node_pt(6);
4733 
4734  brick_face_node_pt[6]=brick_el1_pt->node_pt(9);
4735  brick_face_node_pt[7]=brick_el3_pt->node_pt(1);
4736  brick_face_node_pt[8]=brick_el3_pt->node_pt(3);
4737 
4738  brick_face_node_pt[9]=brick_el0_pt->node_pt(9);
4739  brick_face_node_pt[10]=brick_el0_pt->node_pt(1);
4740  brick_face_node_pt[11]=brick_el1_pt->node_pt(3);
4741 
4742  brick_face_node_pt[12]=brick_el1_pt->node_pt(24);
4743 
4744  brick_face_node_pt[13]=brick_el0_pt->node_pt(10);
4745  brick_face_node_pt[14]=brick_el1_pt->node_pt(21);
4746 
4747  brick_face_node_pt[15]=brick_el1_pt->node_pt(12);
4748  brick_face_node_pt[16]=brick_el3_pt->node_pt(7);
4749 
4750  brick_face_node_pt[17]=brick_el3_pt->node_pt(4);
4751  brick_face_node_pt[18]=brick_el1_pt->node_pt(15);
4752  break;
4753 
4754  }
4755 
4756  // Provide possibility for translation -- may need to add
4757  // face index to this; see ThinLayerBrickOnTetMesh.
4758  Vector<unsigned> translate(19);
4759 
4760  // Initialise with identity mapping
4761  for (unsigned i=0;i<19;i++)
4762  {
4763  translate[i]=i;
4764  }
4765 
4766  // Visit all the nodes on that face
4767  for (unsigned j=0;j<19;j++)
4768  {
4769  // Which node is it?
4770  Node* brick_node_pt=brick_face_node_pt[translate[j]];
4771 
4772  // Get coordinates etc of point from face
4773  Vector<double> s=s_face[j];
4774  Vector<double> zeta(2);
4775  Vector<double> x(3);
4776  face_el_pt->interpolated_zeta(s,zeta);
4777  face_el_pt->interpolated_x(s,x);
4778 
4779 #ifdef PARANOID
4780  // Check that the coordinates match (within tolerance)
4781  double dist=sqrt(pow(brick_node_pt->x(0)-x[0],2)+
4782  pow(brick_node_pt->x(1)-x[1],2)+
4783  pow(brick_node_pt->x(2)-x[2],2));
4785  {
4786  std::ofstream brick0;
4787  std::ofstream brick1;
4788  std::ofstream brick2;
4789  std::ofstream brick3;
4790  brick0.open("full_brick0.dat");
4791  brick1.open("full_brick1.dat");
4792  brick2.open("full_brick2.dat");
4793  brick3.open("full_brick3.dat");
4794  for (unsigned j=0;j<27;j++)
4795  {
4796  brick0 << brick_el0_pt->node_pt(j)->x(0) << " "
4797  << brick_el0_pt->node_pt(j)->x(1) << " "
4798  << brick_el0_pt->node_pt(j)->x(2) << "\n";
4799 
4800  brick1 << brick_el1_pt->node_pt(j)->x(0) << " "
4801  << brick_el1_pt->node_pt(j)->x(1) << " "
4802  << brick_el1_pt->node_pt(j)->x(2) << "\n";
4803 
4804  brick2 << brick_el2_pt->node_pt(j)->x(0) << " "
4805  << brick_el2_pt->node_pt(j)->x(1) << " "
4806  << brick_el2_pt->node_pt(j)->x(2) << "\n";
4807 
4808  brick3 << brick_el3_pt->node_pt(j)->x(0) << " "
4809  << brick_el3_pt->node_pt(j)->x(1) << " "
4810  << brick_el3_pt->node_pt(j)->x(2) << "\n";
4811  }
4812  brick0.close();
4813  brick1.close();
4814  brick2.close();
4815  brick3.close();
4816 
4817  std::ofstream full_face;
4818  full_face.open("full_face.dat");
4819  for (unsigned j=0;j<6;j++)
4820  {
4821  full_face << face_el_pt->node_pt(j)->x(0) << " "
4822  << face_el_pt->node_pt(j)->x(1) << " "
4823  << face_el_pt->node_pt(j)->x(2) << "\n";
4824  }
4825  full_face.close();
4826 
4827  // Get normal sign
4828  int normal_sign=face_el_pt->normal_sign();
4829 
4830  std::ostringstream error_stream;
4831  error_stream
4832  << "During assignment of boundary cordinates, the distance\n"
4833  << "between brick node and reference point in \n"
4834  << "triangular FaceElement is " << dist << " which \n"
4835  << "is bigger than the tolerance defined in \n"
4836  << "BrickFromTetMeshHelper::Face_position_tolerance="
4838  << "If this is tolerable, increase the tolerance \n"
4839  << "(it's defined in a namespace and therefore publically\n"
4840  << "accessible). If not, the Face may be inverted in which \n"
4841  << "case you should re-implement the translation scheme,\n"
4842  << "following the pattern used in the ThinLayerBrickOnTetMesh."
4843  << "\nThe required code fragements are already here but \n"
4844  << "the translation is the unit map.\n"
4845  << "To aid the diagnostics, the files full_brick[0-3].dat\n"
4846  << "contain the coordinates of the 27 nodes in the four\n"
4847  << "bricks associated with the current tet and full_face.dat\n"
4848  << "contains the coordinates of the 6 nodes in the FaceElement"
4849  << "\nfrom which the boundary coordinates are extracted.\n"
4850  << "FYI: The normal_sign of the face is: " << normal_sign
4851  << std::endl;
4852  throw OomphLibError(
4853  error_stream.str(),
4854  OOMPH_CURRENT_FUNCTION,
4855  OOMPH_EXCEPTION_LOCATION);
4856  }
4857 #endif
4858 
4859  // Set boundary stuff
4860  add_boundary_node(b,brick_node_pt);
4861  brick_node_pt->set_coordinates_on_boundary(b,zeta);
4862  }
4863 
4864  // Add appropriate brick elements to boundary lookup scheme
4865  switch (face_index)
4866  {
4867  case 0:
4868  Boundary_element_pt[b].push_back(brick_el1_pt);
4869  Face_index_at_boundary[b].push_back(-2);
4870  Boundary_element_pt[b].push_back(brick_el2_pt);
4871  Face_index_at_boundary[b].push_back(-1);
4872  Boundary_element_pt[b].push_back(brick_el3_pt);
4873  Face_index_at_boundary[b].push_back(-2);
4874  break;
4875 
4876  case 1:
4877  Boundary_element_pt[b].push_back(brick_el0_pt);
4878  Face_index_at_boundary[b].push_back(-1);
4879  Boundary_element_pt[b].push_back(brick_el2_pt);
4880  Face_index_at_boundary[b].push_back(-2);
4881  Boundary_element_pt[b].push_back(brick_el3_pt);
4882  Face_index_at_boundary[b].push_back(-1);
4883  break;
4884 
4885  case 2:
4886  Boundary_element_pt[b].push_back(brick_el0_pt);
4887  Face_index_at_boundary[b].push_back(-3);
4888  Boundary_element_pt[b].push_back(brick_el1_pt);
4889  Face_index_at_boundary[b].push_back(-3);
4890  Boundary_element_pt[b].push_back(brick_el2_pt);
4891  Face_index_at_boundary[b].push_back(-3);
4892  break;
4893 
4894  case 3:
4895  Boundary_element_pt[b].push_back(brick_el0_pt);
4896  Face_index_at_boundary[b].push_back(-2);
4897  Boundary_element_pt[b].push_back(brick_el1_pt);
4898  Face_index_at_boundary[b].push_back(-1);
4899  Boundary_element_pt[b].push_back(brick_el3_pt);
4900  Face_index_at_boundary[b].push_back(-3);
4901  break;
4902  }
4903  // Cleanup
4904  delete face_el_pt;
4905  }
4906  }
4907  // Cleanup
4908  delete face_pt;
4909  }
4910  }
4911 
4912  //Lookup scheme has now been setup
4913  Lookup_for_elements_next_boundary_is_setup=true;
4914 
4915  // Get number of distinct boundaries specified
4916  // in the original xda enumeration.
4917  unsigned n_xda_boundaries=tet_mesh_pt->nxda_boundary();
4918 
4919  // Copy collective IDs across
4920  Boundary_id.resize(n_xda_boundaries);
4921  for (unsigned xda_b=0;xda_b<n_xda_boundaries;xda_b++)
4922  {
4923  Boundary_id[xda_b]=tet_mesh_pt->oomph_lib_boundary_ids(xda_b);
4924  }
4925 
4926 
4927  // Cleanup
4928  for (unsigned e=0;e<4;e++)
4929  {
4930  for (unsigned j=0;j<8;j++)
4931  {
4932  delete dummy_q_el_pt[e]->node_pt(j);
4933  }
4934  delete dummy_q_el_pt[e];
4935  }
4936 
4937  }
4938 
4939  //=======================================================================
4940  /// Build fct: Pass pointer to existing tet mesh and timestepper
4941  /// Specialisation for TetgenMesh<TElement<3,3> >
4942  //=======================================================================
4943  template<class ELEMENT>
4945  TetgenMesh<TElement<3,3> >* tet_mesh_pt,
4946  TimeStepper* time_stepper_pt)
4947  {
4948  // Mesh can only be built with 3D Qelements.
4949  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3, 3);
4950 
4951  // Figure out if the tet mesh is a solid mesh
4952  bool tet_mesh_is_solid_mesh=false;
4953  if (dynamic_cast<SolidFiniteElement*>(tet_mesh_pt->element_pt(0))!=0)
4954  {
4955  tet_mesh_is_solid_mesh=true;
4956  }
4957 
4958  // Setup lookup scheme for local coordinates on triangular faces.
4959  // The local coordinates identify the points on the triangular
4960  // FaceElements on which we place the bottom layer of the
4961  // brick nodes.
4962  Vector<Vector<double> > s_face(19);
4963  for (unsigned i=0;i<19;i++)
4964  {
4965  s_face[i].resize(2);
4966 
4967  switch (i)
4968  {
4969 
4970  // Vertex nodes
4971 
4972  case 0:
4973  s_face[i][0]=1.0;
4974  s_face[i][1]=0.0;
4975  break;
4976 
4977  case 1:
4978  s_face[i][0]=0.0;
4979  s_face[i][1]=1.0;
4980  break;
4981 
4982  case 2:
4983  s_face[i][0]=0.0;
4984  s_face[i][1]=0.0;
4985  break;
4986 
4987  // Midside nodes
4988 
4989  case 3:
4990  s_face[i][0]=0.5;
4991  s_face[i][1]=0.5;
4992  break;
4993 
4994  case 4:
4995  s_face[i][0]=0.0;
4996  s_face[i][1]=0.5;
4997  break;
4998 
4999  case 5:
5000  s_face[i][0]=0.5;
5001  s_face[i][1]=0.0;
5002  break;
5003 
5004 
5005  // Quarter side nodes
5006 
5007  case 6:
5008  s_face[i][0]=0.75;
5009  s_face[i][1]=0.25;
5010  break;
5011 
5012  case 7:
5013  s_face[i][0]=0.25;
5014  s_face[i][1]=0.75;
5015  break;
5016 
5017  case 8:
5018  s_face[i][0]=0.0;
5019  s_face[i][1]=0.75;
5020  break;
5021 
5022  case 9:
5023  s_face[i][0]=0.0;
5024  s_face[i][1]=0.25;
5025  break;
5026 
5027  case 10:
5028  s_face[i][0]=0.25;
5029  s_face[i][1]=0.0;
5030  break;
5031 
5032  case 11:
5033  s_face[i][0]=0.75;
5034  s_face[i][1]=0.0;
5035  break;
5036 
5037  // Central node
5038 
5039  case 12:
5040  s_face[i][0]=1.0/3.0;
5041  s_face[i][1]=1.0/3.0;
5042  break;
5043 
5044 
5045  // Vertical internal midside nodes connecting 2 and 3
5046 
5047  case 13:
5048  s_face[i][0]=5.0/24.0;
5049  s_face[i][1]=5.0/24.0;
5050  break;
5051 
5052  case 14:
5053  s_face[i][0]=5.0/12.0;
5054  s_face[i][1]=5.0/12.0;
5055  break;
5056 
5057  // Internal midside nodes connecting nodes 0 and 4
5058 
5059  case 15:
5060  s_face[i][1]=5.0/24.0;
5061  s_face[i][0]=7.0/12.0; // 1.0-2.0*5.0/24.0;
5062  break;
5063 
5064  case 16:
5065  s_face[i][1]=5.0/12.0;
5066  s_face[i][0]=1.0/6.0; // 1.0-2.0*5.0/12.0;
5067  break;
5068 
5069 
5070  // Internal midside nodes connecting nodes 1 and 5
5071 
5072  case 17:
5073  s_face[i][0]=5.0/24.0;
5074  s_face[i][1]=7.0/12.0; // 1.0-2.0*5.0/24.0;
5075  break;
5076 
5077  case 18:
5078  s_face[i][0]=5.0/12.0;
5079  s_face[i][1]=1.0/6.0; //1.0-2.0*5.0/12.0;
5080  break;
5081 
5082  }
5083  }
5084 
5085  // Set number of boundaries
5086  unsigned nb=tet_mesh_pt->nboundary();
5087  set_nboundary(nb);
5088 
5089  // Get ready for boundary lookup scheme
5090  Boundary_element_pt.resize(nb);
5091  Face_index_at_boundary.resize(nb);
5092 
5093  // Maps to check which nodes have already been done
5094 
5095  // Map that stores the new brick node corresponding to an existing tet node
5096  std::map<Node*,Node*> tet_node_node_pt;
5097 
5098  // Map that stores node on an edge between two brick nodes
5099  std::map<Edge,Node*> brick_edge_node_pt;
5100 
5101  // Map that stores node on face spanned by three tet nodes
5102  std::map<TFace,Node*> tet_face_node_pt;
5103 
5104  // Create the four Dummy bricks:
5105  //------------------------------
5106  Vector<DummyBrickElement*> dummy_q_el_pt(4);
5107  for (unsigned e=0;e<4;e++)
5108  {
5109  dummy_q_el_pt[e]=new DummyBrickElement;
5110  for (unsigned j=0;j<8;j++)
5111  {
5112  dummy_q_el_pt[e]->construct_node(j);
5113  }
5114  }
5115 
5116  // Loop over the elements in the tet mesh
5117  unsigned n_el_tet=tet_mesh_pt->nelement();
5118  for (unsigned e_tet=0; e_tet<n_el_tet; e_tet++)
5119  {
5120  // Cast to ten-noded tet
5121  TElement<3,3>* tet_el_pt=dynamic_cast<TElement<3,3>*>(
5122  tet_mesh_pt->element_pt(e_tet));
5123 
5124 #ifdef PARANOID
5125  if (tet_el_pt==0)
5126  {
5127  std::ostringstream error_stream;
5128  error_stream
5129  << "BrickFromTetMesh can only built from tet mesh containing\n"
5130  << "ten-noded tets.\n";
5131  throw OomphLibError(
5132  error_stream.str(),
5133  OOMPH_CURRENT_FUNCTION,
5134  OOMPH_EXCEPTION_LOCATION);
5135  }
5136 #endif
5137 
5138  // Storage for the centroid node for this tet
5139  Node* centroid_node_pt=0;
5140 
5141  // Internal mid brick-face nodes
5142  Node* top_mid_face_node0_pt=0;
5143  Node* right_mid_face_node0_pt=0;
5144  Node* back_mid_face_node0_pt=0;
5145 
5146  Node* top_mid_face_node1_pt=0;
5147  Node* right_mid_face_node1_pt=0;
5148 
5149  Node* top_mid_face_node2_pt=0;
5150 
5151  // Newly created brick elements
5152  FiniteElement* brick_el0_pt=0;
5153  FiniteElement* brick_el1_pt=0;
5154  FiniteElement* brick_el2_pt=0;
5155  FiniteElement* brick_el3_pt=0;
5156 
5157 
5158  // First brick element is centred at node 0 of tet:
5159  //-------------------------------------------------
5160  {
5161 
5162  // Assign coordinates of dummy element
5163  for (unsigned j=0;j<8;j++)
5164  {
5165  Node* nod_pt=dummy_q_el_pt[0]->node_pt(j);
5166  Vector<double> s_tet(3);
5167  Vector<double> x_tet(3);
5168  switch (j)
5169  {
5170  case 0:
5171  tet_el_pt->local_coordinate_of_node(0,s_tet);
5172  nod_pt->set_value(0,s_tet[0]);
5173  nod_pt->set_value(1,s_tet[1]);
5174  nod_pt->set_value(2,s_tet[2]);
5175  tet_el_pt->interpolated_x(s_tet,x_tet);
5176  nod_pt->x(0)=x_tet[0];
5177  nod_pt->x(1)=x_tet[1];
5178  nod_pt->x(2)=x_tet[2];
5179  break;
5180  case 1:
5181  tet_el_pt->local_coordinate_of_node(4,s_tet);
5182  nod_pt->set_value(0,s_tet[0]);
5183  nod_pt->set_value(1,s_tet[1]);
5184  nod_pt->set_value(2,s_tet[2]);
5185  tet_el_pt->interpolated_x(s_tet,x_tet);
5186  nod_pt->x(0)=x_tet[0];
5187  nod_pt->x(1)=x_tet[1];
5188  nod_pt->x(2)=x_tet[2];
5189  break;
5190  case 2:
5191  tet_el_pt->local_coordinate_of_node(6,s_tet);
5192  nod_pt->set_value(0,s_tet[0]);
5193  nod_pt->set_value(1,s_tet[1]);
5194  nod_pt->set_value(2,s_tet[2]);
5195  tet_el_pt->interpolated_x(s_tet,x_tet);
5196  nod_pt->x(0)=x_tet[0];
5197  nod_pt->x(1)=x_tet[1];
5198  nod_pt->x(2)=x_tet[2];
5199  break;
5200  case 3:
5201  // label 13 in initial sketch: Mid face node on face spanned by
5202  // tet nodes 0,1,3
5203  s_tet[0]=1.0/3.0;
5204  s_tet[1]=1.0/3.0;
5205  s_tet[2]=0.0;
5206  nod_pt->set_value(0,s_tet[0]);
5207  nod_pt->set_value(1,s_tet[1]);
5208  nod_pt->set_value(2,s_tet[2]);
5209  tet_el_pt->interpolated_x(s_tet,x_tet);
5210  nod_pt->x(0)=x_tet[0];
5211  nod_pt->x(1)=x_tet[1];
5212  nod_pt->x(2)=x_tet[2];
5213  break;
5214  case 4:
5215  tet_el_pt->local_coordinate_of_node(5,s_tet);
5216  nod_pt->set_value(0,s_tet[0]);
5217  nod_pt->set_value(1,s_tet[1]);
5218  nod_pt->set_value(2,s_tet[2]);
5219  tet_el_pt->interpolated_x(s_tet,x_tet);
5220  nod_pt->x(0)=x_tet[0];
5221  nod_pt->x(1)=x_tet[1];
5222  nod_pt->x(2)=x_tet[2];
5223  break;
5224  case 5:
5225  // label 11 in initial sketch: Mid face node on face spanned
5226  // by tet nodes 0,1,2
5227  s_tet[0]=1.0/3.0;
5228  s_tet[1]=1.0/3.0;
5229  s_tet[2]=1.0/3.0;
5230  nod_pt->set_value(0,s_tet[0]);
5231  nod_pt->set_value(1,s_tet[1]);
5232  nod_pt->set_value(2,s_tet[2]);
5233  tet_el_pt->interpolated_x(s_tet,x_tet);
5234  nod_pt->x(0)=x_tet[0];
5235  nod_pt->x(1)=x_tet[1];
5236  nod_pt->x(2)=x_tet[2];
5237  break;
5238  case 6:
5239  // label 12 in initial sketch: Mid face node on face
5240  // spanned by tet nodes 0,2,3
5241  s_tet[0]=1.0/3.0;
5242  s_tet[1]=0.0;
5243  s_tet[2]=1.0/3.0;
5244  nod_pt->set_value(0,s_tet[0]);
5245  nod_pt->set_value(1,s_tet[1]);
5246  nod_pt->set_value(2,s_tet[2]);
5247  tet_el_pt->interpolated_x(s_tet,x_tet);
5248  nod_pt->x(0)=x_tet[0];
5249  nod_pt->x(1)=x_tet[1];
5250  nod_pt->x(2)=x_tet[2];
5251  break;
5252  case 7:
5253  // label 14 in initial sketch: Centroid
5254  s_tet[0]=0.25;
5255  s_tet[1]=0.25;
5256  s_tet[2]=0.25;
5257  nod_pt->set_value(0,s_tet[0]);
5258  nod_pt->set_value(1,s_tet[1]);
5259  nod_pt->set_value(2,s_tet[2]);
5260  tet_el_pt->interpolated_x(s_tet,x_tet);
5261  nod_pt->x(0)=x_tet[0];
5262  nod_pt->x(1)=x_tet[1];
5263  nod_pt->x(2)=x_tet[2];
5264  break;
5265  }
5266  }
5267 
5268 
5269  // Create actual zeroth brick element
5270  FiniteElement* el_pt=new ELEMENT;
5271  brick_el0_pt=el_pt;
5272  Element_pt.push_back(el_pt);
5273 
5274  TFace face0(tet_el_pt->node_pt(0),
5275  tet_el_pt->node_pt(1),
5276  tet_el_pt->node_pt(2));
5277 
5278  TFace face1(tet_el_pt->node_pt(0),
5279  tet_el_pt->node_pt(2),
5280  tet_el_pt->node_pt(3));
5281 
5282  TFace face2(tet_el_pt->node_pt(0),
5283  tet_el_pt->node_pt(1),
5284  tet_el_pt->node_pt(3));
5285 
5286 
5287  // Tet vertex nodes along edges emanating from node 0 in brick
5288  Vector<Vector<unsigned> > tet_edge_node(3);
5289  tet_edge_node[0].resize(2);
5290  tet_edge_node[0][0]=4;
5291  tet_edge_node[0][1]=1;
5292  tet_edge_node[1].resize(2);
5293  tet_edge_node[1][0]=6;
5294  tet_edge_node[1][1]=3;
5295  tet_edge_node[2].resize(2);
5296  tet_edge_node[2][0]=5;
5297  tet_edge_node[2][1]=2;
5298 
5299  // Node number of tet vertex that node 0 in brick is centred on
5300  unsigned central_tet_vertex=0;
5301 
5302  Node* tet_node_pt=0;
5303  Node* old_node_pt=0;
5304 
5305  // Corner node
5306  {
5307  unsigned j=0;
5308 
5309  // Need new node?
5310  tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
5311  old_node_pt=tet_node_node_pt[tet_node_pt];
5312  if (old_node_pt==0)
5313  {
5314  Node* new_node_pt=0;
5315  if (tet_node_pt->is_on_boundary())
5316  {
5317  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5318  }
5319  else
5320  {
5321  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5322  }
5323  tet_node_node_pt[tet_node_pt]=new_node_pt;
5324  Node_pt.push_back(new_node_pt);
5325  Vector<double> s(3);
5326  Vector<double> s_tet(3);
5327  Vector<double> x_tet(3);
5328  el_pt->local_coordinate_of_node(j,s);
5329  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5330  tet_el_pt->interpolated_x(s_tet,x_tet);
5331  new_node_pt->x(0)=x_tet[0];
5332  new_node_pt->x(1)=x_tet[1];
5333  new_node_pt->x(2)=x_tet[2];
5334  }
5335  // Node already exists
5336  else
5337  {
5338  el_pt->node_pt(j)=old_node_pt;
5339  }
5340  }
5341 
5342 
5343  // Brick vertex node coindides with mid-edge node on tet edge 0
5344  {
5345  unsigned j=2;
5346 
5347  // Need new node?
5348  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
5349  old_node_pt=tet_node_node_pt[tet_node_pt];
5350  if (old_node_pt==0)
5351  {
5352  Node* new_node_pt=0;
5353  if (tet_node_pt->is_on_boundary())
5354  {
5355  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5356  }
5357  else
5358  {
5359  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5360  }
5361  tet_node_node_pt[tet_node_pt]=new_node_pt;
5362  Node_pt.push_back(new_node_pt);
5363  Vector<double> s(3);
5364  Vector<double> s_tet(3);
5365  Vector<double> x_tet(3);
5366  el_pt->local_coordinate_of_node(j,s);
5367  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5368  tet_el_pt->interpolated_x(s_tet,x_tet);
5369  new_node_pt->x(0)=x_tet[0];
5370  new_node_pt->x(1)=x_tet[1];
5371  new_node_pt->x(2)=x_tet[2];
5372  }
5373  // Node already exists
5374  else
5375  {
5376  el_pt->node_pt(j)=old_node_pt;
5377  }
5378  }
5379 
5380 
5381  // Brick vertex node coindides with mid vertex node of tet edge 1
5382  {
5383  unsigned j=6;
5384 
5385  // Need new node?
5386  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
5387  old_node_pt=tet_node_node_pt[tet_node_pt];
5388  if (old_node_pt==0)
5389  {
5390  Node* new_node_pt=0;
5391  if (tet_node_pt->is_on_boundary())
5392  {
5393  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5394  }
5395  else
5396  {
5397  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5398  }
5399  tet_node_node_pt[tet_node_pt]=new_node_pt;
5400  Node_pt.push_back(new_node_pt);
5401  Vector<double> s(3);
5402  Vector<double> s_tet(3);
5403  Vector<double> x_tet(3);
5404  el_pt->local_coordinate_of_node(j,s);
5405  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5406  tet_el_pt->interpolated_x(s_tet,x_tet);
5407  new_node_pt->x(0)=x_tet[0];
5408  new_node_pt->x(1)=x_tet[1];
5409  new_node_pt->x(2)=x_tet[2];
5410  }
5411  // Node already exists
5412  else
5413  {
5414  el_pt->node_pt(j)=old_node_pt;
5415  }
5416  }
5417 
5418 
5419  // Brick vertex node coindides with mid-vertex node of tet edge 2
5420  {
5421  unsigned j=18;
5422 
5423  // Need new node?
5424  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
5425  old_node_pt=tet_node_node_pt[tet_node_pt];
5426  if (old_node_pt==0)
5427  {
5428  Node* new_node_pt=0;
5429  if (tet_node_pt->is_on_boundary())
5430  {
5431  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5432  }
5433  else
5434  {
5435  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5436  }
5437  tet_node_node_pt[tet_node_pt]=new_node_pt;
5438  Node_pt.push_back(new_node_pt);
5439  Vector<double> s(3);
5440  Vector<double> s_tet(3);
5441  Vector<double> x_tet(3);
5442  el_pt->local_coordinate_of_node(j,s);
5443  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5444  tet_el_pt->interpolated_x(s_tet,x_tet);
5445  new_node_pt->x(0)=x_tet[0];
5446  new_node_pt->x(1)=x_tet[1];
5447  new_node_pt->x(2)=x_tet[2];
5448  }
5449  // Node already exists
5450  else
5451  {
5452  el_pt->node_pt(j)=old_node_pt;
5453  }
5454  }
5455 
5456 
5457 
5458  // Brick vertex node in the middle of tet face0, spanned by
5459  // tet vertices 0, 1, 2. Enumerated "11" in initial sketch.
5460  {
5461  unsigned j=20;
5462 
5463  // Need new node?
5464  old_node_pt=tet_face_node_pt[face0];
5465  if (old_node_pt==0)
5466  {
5467  Node* new_node_pt=0;
5468  if (face0.is_boundary_face())
5469  {
5470  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5471  }
5472  else
5473  {
5474  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5475  }
5476  tet_face_node_pt[face0]=new_node_pt;
5477  Node_pt.push_back(new_node_pt);
5478  Vector<double> s(3);
5479  Vector<double> s_tet(3);
5480  Vector<double> x_tet(3);
5481  el_pt->local_coordinate_of_node(j,s);
5482  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5483  tet_el_pt->interpolated_x(s_tet,x_tet);
5484  new_node_pt->x(0)=x_tet[0];
5485  new_node_pt->x(1)=x_tet[1];
5486  new_node_pt->x(2)=x_tet[2];
5487  }
5488  // Node already exists
5489  else
5490  {
5491  el_pt->node_pt(j)=old_node_pt;
5492  }
5493  }
5494 
5495  // Brick vertex node in the middle of tet face1, spanned by
5496  // tet vertices 0, 2, 3. Enumerated "12" in initial sketch.
5497  {
5498  unsigned j=24;
5499 
5500  // Need new node?
5501  old_node_pt=tet_face_node_pt[face1];
5502  if (old_node_pt==0)
5503  {
5504  Node* new_node_pt=0;
5505  if (face1.is_boundary_face())
5506  {
5507  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5508  }
5509  else
5510  {
5511  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5512  }
5513  tet_face_node_pt[face1]=new_node_pt;
5514  Node_pt.push_back(new_node_pt);
5515  Vector<double> s(3);
5516  Vector<double> s_tet(3);
5517  Vector<double> x_tet(3);
5518  el_pt->local_coordinate_of_node(j,s);
5519  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5520  tet_el_pt->interpolated_x(s_tet,x_tet);
5521  new_node_pt->x(0)=x_tet[0];
5522  new_node_pt->x(1)=x_tet[1];
5523  new_node_pt->x(2)=x_tet[2];
5524  }
5525  // Node already exists
5526  else
5527  {
5528  el_pt->node_pt(j)=old_node_pt;
5529  }
5530  }
5531 
5532  // Brick vertex node in the middle of tet face2, spanned by
5533  // tet vertices 0, 1, 3. Enumerated "13" in initial sketch.
5534  {
5535  unsigned j=8;
5536 
5537  // Need new node?
5538  old_node_pt=tet_face_node_pt[face2];
5539  if (old_node_pt==0)
5540  {
5541  Node* new_node_pt=0;
5542  if (face2.is_boundary_face())
5543  {
5544  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5545  }
5546  else
5547  {
5548  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5549  }
5550  tet_face_node_pt[face2]=new_node_pt;
5551  Node_pt.push_back(new_node_pt);
5552  Vector<double> s(3);
5553  Vector<double> s_tet(3);
5554  Vector<double> x_tet(3);
5555  el_pt->local_coordinate_of_node(j,s);
5556  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5557  tet_el_pt->interpolated_x(s_tet,x_tet);
5558  new_node_pt->x(0)=x_tet[0];
5559  new_node_pt->x(1)=x_tet[1];
5560  new_node_pt->x(2)=x_tet[2];
5561  }
5562  // Node already exists
5563  else
5564  {
5565  el_pt->node_pt(j)=old_node_pt;
5566  }
5567  }
5568 
5569  // Brick vertex node in centroid of tet. Only built for first element.
5570  // Enumerated "13" in initial sketch.
5571  {
5572  unsigned j=26;
5573 
5574  // Always new
5575  {
5576  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5577  centroid_node_pt=new_node_pt;
5578  Node_pt.push_back(new_node_pt);
5579  Vector<double> s(3);
5580  Vector<double> s_tet(3);
5581  Vector<double> x_tet(3);
5582  el_pt->local_coordinate_of_node(j,s);
5583  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5584  tet_el_pt->interpolated_x(s_tet,x_tet);
5585  new_node_pt->x(0)=x_tet[0];
5586  new_node_pt->x(1)=x_tet[1];
5587  new_node_pt->x(2)=x_tet[2];
5588  }
5589  }
5590 
5591 
5592  // Internal brick node -- always built
5593  {
5594  unsigned j=13;
5595 
5596  // Always new
5597  {
5598  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5599  Node_pt.push_back(new_node_pt);
5600  Vector<double> s(3);
5601  Vector<double> s_tet(3);
5602  Vector<double> x_tet(3);
5603  el_pt->local_coordinate_of_node(j,s);
5604  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5605  tet_el_pt->interpolated_x(s_tet,x_tet);
5606  new_node_pt->x(0)=x_tet[0];
5607  new_node_pt->x(1)=x_tet[1];
5608  new_node_pt->x(2)=x_tet[2];
5609  }
5610  }
5611 
5612  // Brick edge node between brick nodes 0 and 2
5613  {
5614  unsigned j=1;
5615 
5616  // Need new node?
5617  Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
5618  old_node_pt=brick_edge_node_pt[edge];
5619  if (old_node_pt==0)
5620  {
5621  Node* new_node_pt=0;
5622  if (edge.is_boundary_edge())
5623  {
5624  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5625  }
5626  else
5627  {
5628  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5629  }
5630  brick_edge_node_pt[edge]=new_node_pt;
5631  Node_pt.push_back(new_node_pt);
5632  Vector<double> s(3);
5633  Vector<double> s_tet(3);
5634  Vector<double> x_tet(3);
5635  el_pt->local_coordinate_of_node(j,s);
5636  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5637  tet_el_pt->interpolated_x(s_tet,x_tet);
5638  new_node_pt->x(0)=x_tet[0];
5639  new_node_pt->x(1)=x_tet[1];
5640  new_node_pt->x(2)=x_tet[2];
5641  }
5642  // Node already exists
5643  else
5644  {
5645  el_pt->node_pt(j)=old_node_pt;
5646  }
5647  }
5648 
5649 
5650  // Brick edge node between brick nodes 0 and 6
5651  {
5652  unsigned j=3;
5653 
5654  // Need new node?
5655  Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
5656  old_node_pt=brick_edge_node_pt[edge];
5657  if (old_node_pt==0)
5658  {
5659  Node* new_node_pt=0;
5660  if (edge.is_boundary_edge())
5661  {
5662  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5663  }
5664  else
5665  {
5666  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5667  }
5668  brick_edge_node_pt[edge]=new_node_pt;
5669  Node_pt.push_back(new_node_pt);
5670  Vector<double> s(3);
5671  Vector<double> s_tet(3);
5672  Vector<double> x_tet(3);
5673  el_pt->local_coordinate_of_node(j,s);
5674  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5675  tet_el_pt->interpolated_x(s_tet,x_tet);
5676  new_node_pt->x(0)=x_tet[0];
5677  new_node_pt->x(1)=x_tet[1];
5678  new_node_pt->x(2)=x_tet[2];
5679  }
5680  // Node already exists
5681  else
5682  {
5683  el_pt->node_pt(j)=old_node_pt;
5684  }
5685  }
5686 
5687  // Brick edge node between brick nodes 2 and 8
5688  {
5689  unsigned j=5;
5690 
5691  // Need new node?
5692  Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
5693  old_node_pt=brick_edge_node_pt[edge];
5694  if (old_node_pt==0)
5695  {
5696  Node* new_node_pt=0;
5697  if (edge.is_boundary_edge())
5698  {
5699  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5700  }
5701  else
5702  {
5703  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5704  }
5705  brick_edge_node_pt[edge]=new_node_pt;
5706  Node_pt.push_back(new_node_pt);
5707  Vector<double> s(3);
5708  Vector<double> s_tet(3);
5709  Vector<double> x_tet(3);
5710  el_pt->local_coordinate_of_node(j,s);
5711  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5712  tet_el_pt->interpolated_x(s_tet,x_tet);
5713  new_node_pt->x(0)=x_tet[0];
5714  new_node_pt->x(1)=x_tet[1];
5715  new_node_pt->x(2)=x_tet[2];
5716  }
5717  // Node already exists
5718  else
5719  {
5720  el_pt->node_pt(j)=old_node_pt;
5721  }
5722  }
5723 
5724  // Brick edge node between brick nodes 6 and 8
5725  {
5726  unsigned j=7;
5727 
5728  // Need new node?
5729  Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
5730  old_node_pt=brick_edge_node_pt[edge];
5731  if (old_node_pt==0)
5732  {
5733  Node* new_node_pt=0;
5734  if (edge.is_boundary_edge())
5735  {
5736  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5737  }
5738  else
5739  {
5740  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5741  }
5742  brick_edge_node_pt[edge]=new_node_pt;
5743  Node_pt.push_back(new_node_pt);
5744  Vector<double> s(3);
5745  Vector<double> s_tet(3);
5746  Vector<double> x_tet(3);
5747  el_pt->local_coordinate_of_node(j,s);
5748  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5749  tet_el_pt->interpolated_x(s_tet,x_tet);
5750  new_node_pt->x(0)=x_tet[0];
5751  new_node_pt->x(1)=x_tet[1];
5752  new_node_pt->x(2)=x_tet[2];
5753  }
5754  // Node already exists
5755  else
5756  {
5757  el_pt->node_pt(j)=old_node_pt;
5758  }
5759  }
5760 
5761  // Brick edge node between brick nodes 18 and 20
5762  {
5763  unsigned j=19;
5764 
5765  // Need new node?
5766  Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
5767  old_node_pt=brick_edge_node_pt[edge];
5768  if (old_node_pt==0)
5769  {
5770  Node* new_node_pt=0;
5771  if (edge.is_boundary_edge())
5772  {
5773  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5774  }
5775  else
5776  {
5777  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5778  }
5779  brick_edge_node_pt[edge]=new_node_pt;
5780  Node_pt.push_back(new_node_pt);
5781  Vector<double> s(3);
5782  Vector<double> s_tet(3);
5783  Vector<double> x_tet(3);
5784  el_pt->local_coordinate_of_node(j,s);
5785  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5786  tet_el_pt->interpolated_x(s_tet,x_tet);
5787  new_node_pt->x(0)=x_tet[0];
5788  new_node_pt->x(1)=x_tet[1];
5789  new_node_pt->x(2)=x_tet[2];
5790  }
5791  // Node already exists
5792  else
5793  {
5794  el_pt->node_pt(j)=old_node_pt;
5795  }
5796  }
5797 
5798 
5799  // Brick edge node between brick nodes 18 and 24
5800  {
5801  unsigned j=21;
5802 
5803  // Need new node?
5804  Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
5805  old_node_pt=brick_edge_node_pt[edge];
5806  if (old_node_pt==0)
5807  {
5808  Node* new_node_pt=0;
5809  if (edge.is_boundary_edge())
5810  {
5811  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5812  }
5813  else
5814  {
5815  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5816  }
5817  brick_edge_node_pt[edge]=new_node_pt;
5818  Node_pt.push_back(new_node_pt);
5819  Vector<double> s(3);
5820  Vector<double> s_tet(3);
5821  Vector<double> x_tet(3);
5822  el_pt->local_coordinate_of_node(j,s);
5823  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5824  tet_el_pt->interpolated_x(s_tet,x_tet);
5825  new_node_pt->x(0)=x_tet[0];
5826  new_node_pt->x(1)=x_tet[1];
5827  new_node_pt->x(2)=x_tet[2];
5828  }
5829  // Node already exists
5830  else
5831  {
5832  el_pt->node_pt(j)=old_node_pt;
5833  }
5834  }
5835 
5836  // Brick edge node between brick nodes 20 and 26
5837  {
5838  unsigned j=23;
5839 
5840  // Need new node?
5841  Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
5842  old_node_pt=brick_edge_node_pt[edge];
5843  if (old_node_pt==0)
5844  {
5845  Node* new_node_pt=0;
5846  if (edge.is_boundary_edge())
5847  {
5848  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5849  }
5850  else
5851  {
5852  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5853  }
5854  brick_edge_node_pt[edge]=new_node_pt;
5855  Node_pt.push_back(new_node_pt);
5856  Vector<double> s(3);
5857  Vector<double> s_tet(3);
5858  Vector<double> x_tet(3);
5859  el_pt->local_coordinate_of_node(j,s);
5860  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5861  tet_el_pt->interpolated_x(s_tet,x_tet);
5862  new_node_pt->x(0)=x_tet[0];
5863  new_node_pt->x(1)=x_tet[1];
5864  new_node_pt->x(2)=x_tet[2];
5865  }
5866  // Node already exists
5867  else
5868  {
5869  el_pt->node_pt(j)=old_node_pt;
5870  }
5871  }
5872 
5873 
5874  // Brick edge node between brick nodes 24 and 26
5875  {
5876  unsigned j=25;
5877 
5878  // Need new node?
5879  Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
5880  old_node_pt=brick_edge_node_pt[edge];
5881  if (old_node_pt==0)
5882  {
5883  Node* new_node_pt=0;
5884  if (edge.is_boundary_edge())
5885  {
5886  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5887  }
5888  else
5889  {
5890  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5891  }
5892  brick_edge_node_pt[edge]=new_node_pt;
5893  Node_pt.push_back(new_node_pt);
5894  Vector<double> s(3);
5895  Vector<double> s_tet(3);
5896  Vector<double> x_tet(3);
5897  el_pt->local_coordinate_of_node(j,s);
5898  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5899  tet_el_pt->interpolated_x(s_tet,x_tet);
5900  new_node_pt->x(0)=x_tet[0];
5901  new_node_pt->x(1)=x_tet[1];
5902  new_node_pt->x(2)=x_tet[2];
5903  }
5904  // Node already exists
5905  else
5906  {
5907  el_pt->node_pt(j)=old_node_pt;
5908  }
5909  }
5910 
5911  // Brick edge node between brick nodes 0 and 18
5912  {
5913  unsigned j=9;
5914 
5915  // Need new node?
5916  Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
5917  old_node_pt=brick_edge_node_pt[edge];
5918  if (old_node_pt==0)
5919  {
5920  Node* new_node_pt=0;
5921  if (edge.is_boundary_edge())
5922  {
5923  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5924  }
5925  else
5926  {
5927  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5928  }
5929  brick_edge_node_pt[edge]=new_node_pt;
5930  Node_pt.push_back(new_node_pt);
5931  Vector<double> s(3);
5932  Vector<double> s_tet(3);
5933  Vector<double> x_tet(3);
5934  el_pt->local_coordinate_of_node(j,s);
5935  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5936  tet_el_pt->interpolated_x(s_tet,x_tet);
5937  new_node_pt->x(0)=x_tet[0];
5938  new_node_pt->x(1)=x_tet[1];
5939  new_node_pt->x(2)=x_tet[2];
5940  }
5941  // Node already exists
5942  else
5943  {
5944  el_pt->node_pt(j)=old_node_pt;
5945  }
5946  }
5947 
5948 
5949  // Brick edge node between brick nodes 2 and 20
5950  {
5951  unsigned j=11;
5952 
5953  // Need new node?
5954  Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
5955  old_node_pt=brick_edge_node_pt[edge];
5956  if (old_node_pt==0)
5957  {
5958  Node* new_node_pt=0;
5959  if (edge.is_boundary_edge())
5960  {
5961  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5962  }
5963  else
5964  {
5965  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5966  }
5967  brick_edge_node_pt[edge]=new_node_pt;
5968  Node_pt.push_back(new_node_pt);
5969  Vector<double> s(3);
5970  Vector<double> s_tet(3);
5971  Vector<double> x_tet(3);
5972  el_pt->local_coordinate_of_node(j,s);
5973  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5974  tet_el_pt->interpolated_x(s_tet,x_tet);
5975  new_node_pt->x(0)=x_tet[0];
5976  new_node_pt->x(1)=x_tet[1];
5977  new_node_pt->x(2)=x_tet[2];
5978  }
5979  // Node already exists
5980  else
5981  {
5982  el_pt->node_pt(j)=old_node_pt;
5983  }
5984  }
5985 
5986 
5987  // Brick edge node between brick nodes 6 and 24
5988  {
5989  unsigned j=15;
5990 
5991  // Need new node?
5992  Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
5993  old_node_pt=brick_edge_node_pt[edge];
5994  if (old_node_pt==0)
5995  {
5996  Node* new_node_pt=0;
5997  if (edge.is_boundary_edge())
5998  {
5999  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6000  }
6001  else
6002  {
6003  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6004  }
6005  brick_edge_node_pt[edge]=new_node_pt;
6006  Node_pt.push_back(new_node_pt);
6007  Vector<double> s(3);
6008  Vector<double> s_tet(3);
6009  Vector<double> x_tet(3);
6010  el_pt->local_coordinate_of_node(j,s);
6011  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6012  tet_el_pt->interpolated_x(s_tet,x_tet);
6013  new_node_pt->x(0)=x_tet[0];
6014  new_node_pt->x(1)=x_tet[1];
6015  new_node_pt->x(2)=x_tet[2];
6016  }
6017  // Node already exists
6018  else
6019  {
6020  el_pt->node_pt(j)=old_node_pt;
6021  }
6022  }
6023 
6024 
6025  // Brick edge node between brick nodes 8 and 26
6026  {
6027  unsigned j=17;
6028 
6029  // Need new node?
6030  Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
6031  old_node_pt=brick_edge_node_pt[edge];
6032  if (old_node_pt==0)
6033  {
6034  Node* new_node_pt=0;
6035  if (edge.is_boundary_edge())
6036  {
6037  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6038  }
6039  else
6040  {
6041  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6042  }
6043  brick_edge_node_pt[edge]=new_node_pt;
6044  Node_pt.push_back(new_node_pt);
6045  Vector<double> s(3);
6046  Vector<double> s_tet(3);
6047  Vector<double> x_tet(3);
6048  el_pt->local_coordinate_of_node(j,s);
6049  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6050  tet_el_pt->interpolated_x(s_tet,x_tet);
6051  new_node_pt->x(0)=x_tet[0];
6052  new_node_pt->x(1)=x_tet[1];
6053  new_node_pt->x(2)=x_tet[2];
6054  }
6055  // Node already exists
6056  else
6057  {
6058  el_pt->node_pt(j)=old_node_pt;
6059  }
6060  }
6061 
6062 
6063  // Mid brick-face node associated with face
6064  // spanned by mid-vertex nodes associated with tet edges 0 and 2
6065  {
6066  unsigned j=10;
6067 
6068  // Need new node?
6069  TFace face(tet_el_pt->node_pt( central_tet_vertex),
6070  tet_el_pt->node_pt(tet_edge_node[0][0]),
6071  tet_el_pt->node_pt(tet_edge_node[2][0]));
6072 
6073  old_node_pt=tet_face_node_pt[face];
6074  if (old_node_pt==0)
6075  {
6076  Node* new_node_pt=0;
6077  if (face.is_boundary_face())
6078  {
6079  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6080  }
6081  else
6082  {
6083  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6084  }
6085  tet_face_node_pt[face]=new_node_pt;
6086  Node_pt.push_back(new_node_pt);
6087  Vector<double> s(3);
6088  Vector<double> s_tet(3);
6089  Vector<double> x_tet(3);
6090  el_pt->local_coordinate_of_node(j,s);
6091  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6092  tet_el_pt->interpolated_x(s_tet,x_tet);
6093  new_node_pt->x(0)=x_tet[0];
6094  new_node_pt->x(1)=x_tet[1];
6095  new_node_pt->x(2)=x_tet[2];
6096  }
6097  // Node already exists
6098  else
6099  {
6100  el_pt->node_pt(j)=old_node_pt;
6101  }
6102  }
6103 
6104 
6105 
6106  // Mid brick-face node associated with face
6107  // spanned by mid-vertex nodes associated with tet edges 1 and 2
6108  {
6109  unsigned j=12;
6110 
6111  // Need new node?
6112  TFace face(tet_el_pt->node_pt(central_tet_vertex),
6113  tet_el_pt->node_pt(tet_edge_node[1][0]),
6114  tet_el_pt->node_pt(tet_edge_node[2][0]));
6115 
6116  old_node_pt=tet_face_node_pt[face];
6117  if (old_node_pt==0)
6118  {
6119  Node* new_node_pt=0;
6120  if (face.is_boundary_face())
6121  {
6122  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6123  }
6124  else
6125  {
6126  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6127  }
6128  tet_face_node_pt[face]=new_node_pt;
6129  Node_pt.push_back(new_node_pt);
6130  Vector<double> s(3);
6131  Vector<double> s_tet(3);
6132  Vector<double> x_tet(3);
6133  el_pt->local_coordinate_of_node(j,s);
6134  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6135  tet_el_pt->interpolated_x(s_tet,x_tet);
6136  new_node_pt->x(0)=x_tet[0];
6137  new_node_pt->x(1)=x_tet[1];
6138  new_node_pt->x(2)=x_tet[2];
6139  }
6140  // Node already exists
6141  else
6142  {
6143  el_pt->node_pt(j)=old_node_pt;
6144  }
6145  }
6146 
6147 
6148 
6149  // Mid brick-face node associated with face
6150  // spanned by mid-vertex nodes associated with tet edges 0 and 1
6151  {
6152  unsigned j=4;
6153 
6154  // Need new node?
6155  TFace face(tet_el_pt->node_pt(central_tet_vertex),
6156  tet_el_pt->node_pt(tet_edge_node[0][0]),
6157  tet_el_pt->node_pt(tet_edge_node[1][0]));
6158 
6159  old_node_pt=tet_face_node_pt[face];
6160  if (old_node_pt==0)
6161  {
6162  Node* new_node_pt=0;
6163  if (face.is_boundary_face())
6164  {
6165  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6166  }
6167  else
6168  {
6169  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6170  }
6171  tet_face_node_pt[face]=new_node_pt;
6172  Node_pt.push_back(new_node_pt);
6173  Vector<double> s(3);
6174  Vector<double> s_tet(3);
6175  Vector<double> x_tet(3);
6176  el_pt->local_coordinate_of_node(j,s);
6177  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6178  tet_el_pt->interpolated_x(s_tet,x_tet);
6179  new_node_pt->x(0)=x_tet[0];
6180  new_node_pt->x(1)=x_tet[1];
6181  new_node_pt->x(2)=x_tet[2];
6182  }
6183  // Node already exists
6184  else
6185  {
6186  el_pt->node_pt(j)=old_node_pt;
6187  }
6188  }
6189 
6190 
6191  // Top mid brick-face node -- only built by first element
6192  {
6193  unsigned j=22;
6194  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6195  Node_pt.push_back(new_node_pt);
6196  Vector<double> s(3);
6197  Vector<double> s_tet(3);
6198  Vector<double> x_tet(3);
6199  el_pt->local_coordinate_of_node(j,s);
6200  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6201  top_mid_face_node0_pt=new_node_pt;
6202  tet_el_pt->interpolated_x(s_tet,x_tet);
6203  new_node_pt->x(0)=x_tet[0];
6204  new_node_pt->x(1)=x_tet[1];
6205  new_node_pt->x(2)=x_tet[2];
6206  }
6207 
6208 
6209 
6210  // Right mid brick-face node -- only built by first element
6211  {
6212  unsigned j=14;
6213  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6214  Node_pt.push_back(new_node_pt);
6215  Vector<double> s(3);
6216  Vector<double> s_tet(3);
6217  Vector<double> x_tet(3);
6218  el_pt->local_coordinate_of_node(j,s);
6219  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6220  right_mid_face_node0_pt=new_node_pt;
6221  tet_el_pt->interpolated_x(s_tet,x_tet);
6222  new_node_pt->x(0)=x_tet[0];
6223  new_node_pt->x(1)=x_tet[1];
6224  new_node_pt->x(2)=x_tet[2];
6225  }
6226 
6227 
6228  // Back mid brick-face node -- only built by first element
6229  {
6230  unsigned j=16;
6231  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6232  Node_pt.push_back(new_node_pt);
6233  Vector<double> s(3);
6234  Vector<double> s_tet(3);
6235  Vector<double> x_tet(3);
6236  el_pt->local_coordinate_of_node(j,s);
6237  dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6238  back_mid_face_node0_pt=new_node_pt;
6239  tet_el_pt->interpolated_x(s_tet,x_tet);
6240  new_node_pt->x(0)=x_tet[0];
6241  new_node_pt->x(1)=x_tet[1];
6242  new_node_pt->x(2)=x_tet[2];
6243  }
6244 
6245  }
6246 
6247 
6248  // Second brick element is centred at node 1 of tet:
6249  //--------------------------------------------------
6250  {
6251  // Assign coordinates of dummy element
6252  for (unsigned j=0;j<8;j++)
6253  {
6254  Node* nod_pt=dummy_q_el_pt[1]->node_pt(j);
6255  Vector<double> s_tet(3);
6256  Vector<double> x_tet(3);
6257  switch (j)
6258  {
6259  case 0:
6260  tet_el_pt->local_coordinate_of_node(1,s_tet);
6261  nod_pt->set_value(0,s_tet[0]);
6262  nod_pt->set_value(1,s_tet[1]);
6263  nod_pt->set_value(2,s_tet[2]);
6264  tet_el_pt->interpolated_x(s_tet,x_tet);
6265  nod_pt->x(0)=x_tet[0];
6266  nod_pt->x(1)=x_tet[1];
6267  nod_pt->x(2)=x_tet[2];
6268  break;
6269  case 1:
6270  tet_el_pt->local_coordinate_of_node(9,s_tet);
6271  nod_pt->set_value(0,s_tet[0]);
6272  nod_pt->set_value(1,s_tet[1]);
6273  nod_pt->set_value(2,s_tet[2]);
6274  tet_el_pt->interpolated_x(s_tet,x_tet);
6275  nod_pt->x(0)=x_tet[0];
6276  nod_pt->x(1)=x_tet[1];
6277  nod_pt->x(2)=x_tet[2];
6278  break;
6279  case 2:
6280  tet_el_pt->local_coordinate_of_node(4,s_tet);
6281  nod_pt->set_value(0,s_tet[0]);
6282  nod_pt->set_value(1,s_tet[1]);
6283  nod_pt->set_value(2,s_tet[2]);
6284  tet_el_pt->interpolated_x(s_tet,x_tet);
6285  nod_pt->x(0)=x_tet[0];
6286  nod_pt->x(1)=x_tet[1];
6287  nod_pt->x(2)=x_tet[2];
6288  break;
6289  case 3:
6290  // label 13 in initial sketch: Mid face node on face
6291  // spanned by tet nodes 0,1,3
6292  s_tet[0]=1.0/3.0;
6293  s_tet[1]=1.0/3.0;
6294  s_tet[2]=0.0;
6295  nod_pt->set_value(0,s_tet[0]);
6296  nod_pt->set_value(1,s_tet[1]);
6297  nod_pt->set_value(2,s_tet[2]);
6298  tet_el_pt->interpolated_x(s_tet,x_tet);
6299  nod_pt->x(0)=x_tet[0];
6300  nod_pt->x(1)=x_tet[1];
6301  nod_pt->x(2)=x_tet[2];
6302  break;
6303  case 4:
6304  tet_el_pt->local_coordinate_of_node(7,s_tet);
6305  nod_pt->set_value(0,s_tet[0]);
6306  nod_pt->set_value(1,s_tet[1]);
6307  nod_pt->set_value(2,s_tet[2]);
6308  tet_el_pt->interpolated_x(s_tet,x_tet);
6309  nod_pt->x(0)=x_tet[0];
6310  nod_pt->x(1)=x_tet[1];
6311  nod_pt->x(2)=x_tet[2];
6312  break;
6313  case 5:
6314  // label 10 in initial sketch: Mid face node on face
6315  // spanned by tet nodes 1,2,3
6316  s_tet[0]=0.0;
6317  s_tet[1]=1.0/3.0;
6318  s_tet[2]=1.0/3.0;
6319  nod_pt->set_value(0,s_tet[0]);
6320  nod_pt->set_value(1,s_tet[1]);
6321  nod_pt->set_value(2,s_tet[2]);
6322  tet_el_pt->interpolated_x(s_tet,x_tet);
6323  nod_pt->x(0)=x_tet[0];
6324  nod_pt->x(1)=x_tet[1];
6325  nod_pt->x(2)=x_tet[2];
6326  break;
6327  case 6:
6328  // label 11 in initial sketch: Mid face node on face
6329  // spanned by tet nodes 0,1,2
6330  s_tet[0]=1.0/3.0;
6331  s_tet[1]=1.0/3.0;
6332  s_tet[2]=1.0/3.0;
6333  nod_pt->set_value(0,s_tet[0]);
6334  nod_pt->set_value(1,s_tet[1]);
6335  nod_pt->set_value(2,s_tet[2]);
6336  tet_el_pt->interpolated_x(s_tet,x_tet);
6337  nod_pt->x(0)=x_tet[0];
6338  nod_pt->x(1)=x_tet[1];
6339  nod_pt->x(2)=x_tet[2];
6340  break;
6341  case 7:
6342  // label 14 in initial sketch: Centroid
6343  s_tet[0]=0.25;
6344  s_tet[1]=0.25;
6345  s_tet[2]=0.25;
6346  nod_pt->set_value(0,s_tet[0]);
6347  nod_pt->set_value(1,s_tet[1]);
6348  nod_pt->set_value(2,s_tet[2]);
6349  tet_el_pt->interpolated_x(s_tet,x_tet);
6350  nod_pt->x(0)=x_tet[0];
6351  nod_pt->x(1)=x_tet[1];
6352  nod_pt->x(2)=x_tet[2];
6353  break;
6354  }
6355  }
6356 
6357 
6358  // Create actual first brick element
6359  FiniteElement* el_pt=new ELEMENT;
6360  brick_el1_pt=el_pt;
6361  Element_pt.push_back(el_pt);
6362 
6363  TFace face0(tet_el_pt->node_pt(1),
6364  tet_el_pt->node_pt(3),
6365  tet_el_pt->node_pt(2));
6366 
6367  TFace face1(tet_el_pt->node_pt(1),
6368  tet_el_pt->node_pt(0),
6369  tet_el_pt->node_pt(2));
6370 
6371  TFace face2(tet_el_pt->node_pt(1),
6372  tet_el_pt->node_pt(0),
6373  tet_el_pt->node_pt(3));
6374 
6375  // Tet vertex nodes along edges emanating from node 0 in brick
6376  Vector<Vector<unsigned> > tet_edge_node(3);
6377  tet_edge_node[0].resize(2);
6378  tet_edge_node[0][0]=9;
6379  tet_edge_node[0][1]=3;
6380  tet_edge_node[1].resize(2);
6381  tet_edge_node[1][0]=4;
6382  tet_edge_node[1][1]=0;
6383  tet_edge_node[2].resize(2);
6384  tet_edge_node[2][0]=7;
6385  tet_edge_node[2][1]=2;
6386 
6387  // Node number of tet vertex that node 0 in brick is centred on
6388  unsigned central_tet_vertex=1;
6389 
6390  Node* tet_node_pt=0;
6391  Node* old_node_pt=0;
6392 
6393  // Corner node
6394  {
6395  unsigned j=0;
6396 
6397  // Need new node?
6398  tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
6399  old_node_pt=tet_node_node_pt[tet_node_pt];
6400  if (old_node_pt==0)
6401  {
6402  Node* new_node_pt=0;
6403  if (tet_node_pt->is_on_boundary())
6404  {
6405  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6406  }
6407  else
6408  {
6409  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6410  }
6411  tet_node_node_pt[tet_node_pt]=new_node_pt;
6412  Node_pt.push_back(new_node_pt);
6413  Vector<double> s(3);
6414  Vector<double> s_tet(3);
6415  Vector<double> x_tet(3);
6416  el_pt->local_coordinate_of_node(j,s);
6417  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6418  tet_el_pt->interpolated_x(s_tet,x_tet);
6419  new_node_pt->x(0)=x_tet[0];
6420  new_node_pt->x(1)=x_tet[1];
6421  new_node_pt->x(2)=x_tet[2];
6422  }
6423  // Node already exists
6424  else
6425  {
6426  el_pt->node_pt(j)=old_node_pt;
6427  }
6428  }
6429 
6430 
6431  // Brick vertex node coindides with mid-edge node on tet edge 0
6432  {
6433  unsigned j=2;
6434 
6435  // Need new node?
6436  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
6437  old_node_pt=tet_node_node_pt[tet_node_pt];
6438  if (old_node_pt==0)
6439  {
6440  Node* new_node_pt=0;
6441  if (tet_node_pt->is_on_boundary())
6442  {
6443  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6444  }
6445  else
6446  {
6447  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6448  }
6449  tet_node_node_pt[tet_node_pt]=new_node_pt;
6450  Node_pt.push_back(new_node_pt);
6451  Vector<double> s(3);
6452  Vector<double> s_tet(3);
6453  Vector<double> x_tet(3);
6454  el_pt->local_coordinate_of_node(j,s);
6455  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6456  tet_el_pt->interpolated_x(s_tet,x_tet);
6457  new_node_pt->x(0)=x_tet[0];
6458  new_node_pt->x(1)=x_tet[1];
6459  new_node_pt->x(2)=x_tet[2];
6460  }
6461  // Node already exists
6462  else
6463  {
6464  el_pt->node_pt(j)=old_node_pt;
6465  }
6466  }
6467 
6468 
6469  // Brick vertex node coindides with mid vertex node of tet edge 1
6470  {
6471  unsigned j=6;
6472 
6473  // Need new node?
6474  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
6475  old_node_pt=tet_node_node_pt[tet_node_pt];
6476  if (old_node_pt==0)
6477  {
6478  Node* new_node_pt=0;
6479  if (tet_node_pt->is_on_boundary())
6480  {
6481  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6482  }
6483  else
6484  {
6485  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6486  }
6487  tet_node_node_pt[tet_node_pt]=new_node_pt;
6488  Node_pt.push_back(new_node_pt);
6489  Vector<double> s(3);
6490  Vector<double> s_tet(3);
6491  Vector<double> x_tet(3);
6492  el_pt->local_coordinate_of_node(j,s);
6493  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6494  tet_el_pt->interpolated_x(s_tet,x_tet);
6495  new_node_pt->x(0)=x_tet[0];
6496  new_node_pt->x(1)=x_tet[1];
6497  new_node_pt->x(2)=x_tet[2];
6498  }
6499  // Node already exists
6500  else
6501  {
6502  el_pt->node_pt(j)=old_node_pt;
6503  }
6504  }
6505 
6506 
6507  // Brick vertex node coindides with mid-vertex node of tet edge 2
6508  {
6509  unsigned j=18;
6510 
6511  // Need new node?
6512  tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
6513  old_node_pt=tet_node_node_pt[tet_node_pt];
6514  if (old_node_pt==0)
6515  {
6516  Node* new_node_pt=0;
6517  if (tet_node_pt->is_on_boundary())
6518  {
6519  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6520  }
6521  else
6522  {
6523  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6524  }
6525  tet_node_node_pt[tet_node_pt]=new_node_pt;
6526  Node_pt.push_back(new_node_pt);
6527  Vector<double> s(3);
6528  Vector<double> s_tet(3);
6529  Vector<double> x_tet(3);
6530  el_pt->local_coordinate_of_node(j,s);
6531  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6532  tet_el_pt->interpolated_x(s_tet,x_tet);
6533  new_node_pt->x(0)=x_tet[0];
6534  new_node_pt->x(1)=x_tet[1];
6535  new_node_pt->x(2)=x_tet[2];
6536  }
6537  // Node already exists
6538  else
6539  {
6540  el_pt->node_pt(j)=old_node_pt;
6541  }
6542  }
6543 
6544 
6545 
6546  // Brick vertex node in the middle of tet face0
6547  {
6548  unsigned j=20;
6549 
6550  // Need new node?
6551  old_node_pt=tet_face_node_pt[face0];
6552  if (old_node_pt==0)
6553  {
6554  Node* new_node_pt=0;
6555  if (face0.is_boundary_face())
6556  {
6557  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6558  }
6559  else
6560  {
6561  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6562  }
6563  tet_face_node_pt[face0]=new_node_pt;
6564  Node_pt.push_back(new_node_pt);
6565  Vector<double> s(3);
6566  Vector<double> s_tet(3);
6567  Vector<double> x_tet(3);
6568  el_pt->local_coordinate_of_node(j,s);
6569  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6570  tet_el_pt->interpolated_x(s_tet,x_tet);
6571  new_node_pt->x(0)=x_tet[0];
6572  new_node_pt->x(1)=x_tet[1];
6573  new_node_pt->x(2)=x_tet[2];
6574  }
6575  // Node already exists
6576  else
6577  {
6578  el_pt->node_pt(j)=old_node_pt;
6579  }
6580  }
6581 
6582  // Brick vertex node in the middle of tet face1
6583  {
6584  unsigned j=24;
6585 
6586  // Need new node?
6587  old_node_pt=tet_face_node_pt[face1];
6588  if (old_node_pt==0)
6589  {
6590  Node* new_node_pt=0;
6591  if (face1.is_boundary_face())
6592  {
6593  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6594  }
6595  else
6596  {
6597  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6598  }
6599  tet_face_node_pt[face1]=new_node_pt;
6600  Node_pt.push_back(new_node_pt);
6601  Vector<double> s(3);
6602  Vector<double> s_tet(3);
6603  Vector<double> x_tet(3);
6604  el_pt->local_coordinate_of_node(j,s);
6605  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6606  tet_el_pt->interpolated_x(s_tet,x_tet);
6607  new_node_pt->x(0)=x_tet[0];
6608  new_node_pt->x(1)=x_tet[1];
6609  new_node_pt->x(2)=x_tet[2];
6610  }
6611  // Node already exists
6612  else
6613  {
6614  el_pt->node_pt(j)=old_node_pt;
6615  }
6616  }
6617 
6618  // Brick vertex node in the middle of face2
6619  {
6620  unsigned j=8;
6621 
6622  // Need new node?
6623  old_node_pt=tet_face_node_pt[face2];
6624  if (old_node_pt==0)
6625  {
6626  Node* new_node_pt=0;
6627  if (face2.is_boundary_face())
6628  {
6629  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6630  }
6631  else
6632  {
6633  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6634  }
6635  tet_face_node_pt[face2]=new_node_pt;
6636  Node_pt.push_back(new_node_pt);
6637  Vector<double> s(3);
6638  Vector<double> s_tet(3);
6639  Vector<double> x_tet(3);
6640  el_pt->local_coordinate_of_node(j,s);
6641  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6642  tet_el_pt->interpolated_x(s_tet,x_tet);
6643  new_node_pt->x(0)=x_tet[0];
6644  new_node_pt->x(1)=x_tet[1];
6645  new_node_pt->x(2)=x_tet[2];
6646  }
6647  // Node already exists
6648  else
6649  {
6650  el_pt->node_pt(j)=old_node_pt;
6651  }
6652  }
6653 
6654 
6655  // Brick vertex node in centroid of tet. Only built for first element.
6656  // Enumerated "13" in initial sketch.
6657  {
6658  unsigned j=26;
6659 
6660  // Always copied
6661  el_pt->node_pt(j)=centroid_node_pt;
6662  }
6663 
6664 
6665  // Internal brick node -- always built
6666  {
6667  unsigned j=13;
6668 
6669  // Always new
6670  {
6671  Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6672  Node_pt.push_back(new_node_pt);
6673  Vector<double> s(3);
6674  Vector<double> s_tet(3);
6675  Vector<double> x_tet(3);
6676  el_pt->local_coordinate_of_node(j,s);
6677  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6678  tet_el_pt->interpolated_x(s_tet,x_tet);
6679  new_node_pt->x(0)=x_tet[0];
6680  new_node_pt->x(1)=x_tet[1];
6681  new_node_pt->x(2)=x_tet[2];
6682  }
6683  }
6684 
6685 
6686  // Brick edge node between brick nodes 0 and 2
6687  {
6688  unsigned j=1;
6689 
6690  // Need new node?
6691  Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
6692  old_node_pt=brick_edge_node_pt[edge];
6693  if (old_node_pt==0)
6694  {
6695  Node* new_node_pt=0;
6696  if (edge.is_boundary_edge())
6697  {
6698  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6699  }
6700  else
6701  {
6702  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6703  }
6704  brick_edge_node_pt[edge]=new_node_pt;
6705  Node_pt.push_back(new_node_pt);
6706  Vector<double> s(3);
6707  Vector<double> s_tet(3);
6708  Vector<double> x_tet(3);
6709  el_pt->local_coordinate_of_node(j,s);
6710  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6711  tet_el_pt->interpolated_x(s_tet,x_tet);
6712  new_node_pt->x(0)=x_tet[0];
6713  new_node_pt->x(1)=x_tet[1];
6714  new_node_pt->x(2)=x_tet[2];
6715  }
6716  // Node already exists
6717  else
6718  {
6719  el_pt->node_pt(j)=old_node_pt;
6720  }
6721  }
6722 
6723 
6724  // Brick edge node between brick nodes 0 and 6
6725  {
6726  unsigned j=3;
6727 
6728  // Need new node?
6729  Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
6730  old_node_pt=brick_edge_node_pt[edge];
6731  if (old_node_pt==0)
6732  {
6733  Node* new_node_pt=0;
6734  if (edge.is_boundary_edge())
6735  {
6736  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6737  }
6738  else
6739  {
6740  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6741  }
6742  brick_edge_node_pt[edge]=new_node_pt;
6743  Node_pt.push_back(new_node_pt);
6744  Vector<double> s(3);
6745  Vector<double> s_tet(3);
6746  Vector<double> x_tet(3);
6747  el_pt->local_coordinate_of_node(j,s);
6748  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6749  tet_el_pt->interpolated_x(s_tet,x_tet);
6750  new_node_pt->x(0)=x_tet[0];
6751  new_node_pt->x(1)=x_tet[1];
6752  new_node_pt->x(2)=x_tet[2];
6753  }
6754  // Node already exists
6755  else
6756  {
6757  el_pt->node_pt(j)=old_node_pt;
6758  }
6759  }
6760 
6761 
6762  // Brick edge node between brick nodes 2 and 8
6763  {
6764  unsigned j=5;
6765 
6766  // Need new node?
6767  Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
6768  old_node_pt=brick_edge_node_pt[edge];
6769  if (old_node_pt==0)
6770  {
6771  Node* new_node_pt=0;
6772  if (edge.is_boundary_edge())
6773  {
6774  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6775  }
6776  else
6777  {
6778  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6779  }
6780  brick_edge_node_pt[edge]=new_node_pt;
6781  Node_pt.push_back(new_node_pt);
6782  Vector<double> s(3);
6783  Vector<double> s_tet(3);
6784  Vector<double> x_tet(3);
6785  el_pt->local_coordinate_of_node(j,s);
6786  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6787  tet_el_pt->interpolated_x(s_tet,x_tet);
6788  new_node_pt->x(0)=x_tet[0];
6789  new_node_pt->x(1)=x_tet[1];
6790  new_node_pt->x(2)=x_tet[2];
6791  }
6792  // Node already exists
6793  else
6794  {
6795  el_pt->node_pt(j)=old_node_pt;
6796  }
6797  }
6798 
6799  // Brick edge node between brick nodes 6 and 8
6800  {
6801  unsigned j=7;
6802 
6803  // Need new node?
6804  Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
6805  old_node_pt=brick_edge_node_pt[edge];
6806  if (old_node_pt==0)
6807  {
6808  Node* new_node_pt=0;
6809  if (edge.is_boundary_edge())
6810  {
6811  new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6812  }
6813  else
6814  {
6815  new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6816  }
6817  brick_edge_node_pt[edge]=new_node_pt;
6818  Node_pt.push_back(new_node_pt);
6819  Vector<double> s(3);
6820  Vector<double> s_tet(3);
6821  Vector<double> x_tet(3);
6822  el_pt->local_coordinate_of_node(j,s);
6823  dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6824  tet_el_pt->interpolated_x(s_tet,x_tet);
6825  new_node_pt->x(0)=x_tet[0];
6826  new_node_pt->x(1)=x_tet[1];
6827  new_node_pt->x(2)=x_tet[2];
6828  }
6829  // Node already exists
6830  else
6831  {
6832  el_pt->node_pt(j)=old_node_pt;
6833  }
6834  }
6835 
6836  // Brick edge node between brick nodes 18 and 20
6837  {
6838