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