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-2022 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
32namespace 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);