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-2025 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_HEADER
27#define OOMPH_BRICK_FROM_TET_MESH_TEMPLATE_HEADER
28
29#ifndef OOMPH_BRICK_FROM_TET_MESH_HEADER
30#error __FILE__ should only be included from brick_from_tet_mesh.h.
31#endif // OOMPH_BRICK_FROM_TET_MESH_HEADER
32
33namespace oomph
34{
35 //=======================================================================
36 /// Build fct: Pass pointer to existing tet mesh and timestepper
37 /// Specialisation for XdaTetMesh<TElement<3,3> >
38 //=======================================================================
39 template<class ELEMENT>
42 {
43 // Mesh can only be built with 3D Qelements.
44 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3, 3);
45
46 // Figure out if the tet mesh is a solid mesh
47 bool tet_mesh_is_solid_mesh = false;
48 if (dynamic_cast<SolidFiniteElement*>(tet_mesh_pt->element_pt(0)) != 0)
49 {
51 }
52
53 // Setup lookup scheme for local coordinates on triangular faces.
54 // The local coordinates identify the points on the triangular
55 // FaceElements on which we place the bottom layer of the
56 // brick nodes.
58 for (unsigned i = 0; i < 19; i++)
59 {
60 s_face[i].resize(2);
61
62 switch (i)
63 {
64 // Vertex nodes
65
66 case 0:
67 s_face[i][0] = 1.0;
68 s_face[i][1] = 0.0;
69 break;
70
71 case 1:
72 s_face[i][0] = 0.0;
73 s_face[i][1] = 1.0;
74 break;
75
76 case 2:
77 s_face[i][0] = 0.0;
78 s_face[i][1] = 0.0;
79 break;
80
81 // Midside nodes
82
83 case 3:
84 s_face[i][0] = 0.5;
85 s_face[i][1] = 0.5;
86 break;
87
88 case 4:
89 s_face[i][0] = 0.0;
90 s_face[i][1] = 0.5;
91 break;
92
93 case 5:
94 s_face[i][0] = 0.5;
95 s_face[i][1] = 0.0;
96 break;
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 // Vertical internal midside nodes connecting 2 and 3
138
139 case 13:
140 s_face[i][0] = 5.0 / 24.0;
141 s_face[i][1] = 5.0 / 24.0;
142 break;
143
144 case 14:
145 s_face[i][0] = 5.0 / 12.0;
146 s_face[i][1] = 5.0 / 12.0;
147 break;
148
149 // Internal midside nodes connecting nodes 0 and 4
150
151 case 15:
152 s_face[i][1] = 5.0 / 24.0;
153 s_face[i][0] = 7.0 / 12.0; // 1.0-2.0*5.0/24.0;
154 break;
155
156 case 16:
157 s_face[i][1] = 5.0 / 12.0;
158 s_face[i][0] = 1.0 / 6.0; // 1.0-2.0*5.0/12.0;
159 break;
160
161 // Internal midside nodes connecting nodes 1 and 5
162
163 case 17:
164 s_face[i][0] = 5.0 / 24.0;
165 s_face[i][1] = 7.0 / 12.0; // 1.0-2.0*5.0/24.0;
166 break;
167
168 case 18:
169 s_face[i][0] = 5.0 / 12.0;
170 s_face[i][1] = 1.0 / 6.0; // 1.0-2.0*5.0/12.0;
171 break;
172 }
173 }
174
175 // Set number of boundaries
176 unsigned nb = tet_mesh_pt->nboundary();
177 set_nboundary(nb);
178
179 // Get ready for boundary lookup scheme
180 Boundary_element_pt.resize(nb);
181 Face_index_at_boundary.resize(nb);
182
183 // Maps to check which nodes have already been done
184
185 // Map that stores the new brick node corresponding to an existing tet node
186 std::map<Node*, Node*> tet_node_node_pt;
187
188 // Map that stores node on an edge between two brick nodes
189 std::map<Edge, Node*> brick_edge_node_pt;
190
191 // Map that stores node on face spanned by three tet nodes
192 std::map<TFace, Node*> tet_face_node_pt;
193
194 // Create the four Dummy bricks:
195 //------------------------------
197 for (unsigned e = 0; e < 4; e++)
198 {
200 for (unsigned j = 0; j < 8; j++)
201 {
203 }
204 }
205
206 // Loop over the elements in the tet mesh
207 unsigned n_el_tet = tet_mesh_pt->nelement();
208 for (unsigned e_tet = 0; e_tet < n_el_tet; e_tet++)
209 {
210 // Cast to ten-noded tet
212 dynamic_cast<TElement<3, 3>*>(tet_mesh_pt->element_pt(e_tet));
213
214#ifdef PARANOID
215 if (tet_el_pt == 0)
216 {
217 std::ostringstream error_stream;
219 << "BrickFromTetMesh can only built from tet mesh containing\n"
220 << "ten-noded tets.\n";
221 throw OomphLibError(
223 }
224#endif
225
226 // Storage for the centroid node for this tet
228
229 // Internal mid brick-face nodes
233
236
238
239 // Newly created brick elements
244
245 // First brick element is centred at node 0 of tet:
246 //-------------------------------------------------
247 {
248 // Assign coordinates of dummy element
249 for (unsigned j = 0; j < 8; j++)
250 {
254 switch (j)
255 {
256 case 0:
258 nod_pt->set_value(0, s_tet[0]);
259 nod_pt->set_value(1, s_tet[1]);
260 nod_pt->set_value(2, s_tet[2]);
262 nod_pt->x(0) = x_tet[0];
263 nod_pt->x(1) = x_tet[1];
264 nod_pt->x(2) = x_tet[2];
265 break;
266 case 1:
268 nod_pt->set_value(0, s_tet[0]);
269 nod_pt->set_value(1, s_tet[1]);
270 nod_pt->set_value(2, s_tet[2]);
272 nod_pt->x(0) = x_tet[0];
273 nod_pt->x(1) = x_tet[1];
274 nod_pt->x(2) = x_tet[2];
275 break;
276 case 2:
278 nod_pt->set_value(0, s_tet[0]);
279 nod_pt->set_value(1, s_tet[1]);
280 nod_pt->set_value(2, s_tet[2]);
282 nod_pt->x(0) = x_tet[0];
283 nod_pt->x(1) = x_tet[1];
284 nod_pt->x(2) = x_tet[2];
285 break;
286 case 3:
287 // label 13 in initial sketch: Mid face node on face spanned by
288 // tet nodes 0,1,3
289 s_tet[0] = 1.0 / 3.0;
290 s_tet[1] = 1.0 / 3.0;
291 s_tet[2] = 0.0;
292 nod_pt->set_value(0, s_tet[0]);
293 nod_pt->set_value(1, s_tet[1]);
294 nod_pt->set_value(2, s_tet[2]);
296 nod_pt->x(0) = x_tet[0];
297 nod_pt->x(1) = x_tet[1];
298 nod_pt->x(2) = x_tet[2];
299 break;
300 case 4:
302 nod_pt->set_value(0, s_tet[0]);
303 nod_pt->set_value(1, s_tet[1]);
304 nod_pt->set_value(2, s_tet[2]);
306 nod_pt->x(0) = x_tet[0];
307 nod_pt->x(1) = x_tet[1];
308 nod_pt->x(2) = x_tet[2];
309 break;
310 case 5:
311 // label 11 in initial sketch: Mid face node on face spanned
312 // by tet nodes 0,1,2
313 s_tet[0] = 1.0 / 3.0;
314 s_tet[1] = 1.0 / 3.0;
315 s_tet[2] = 1.0 / 3.0;
316 nod_pt->set_value(0, s_tet[0]);
317 nod_pt->set_value(1, s_tet[1]);
318 nod_pt->set_value(2, s_tet[2]);
320 nod_pt->x(0) = x_tet[0];
321 nod_pt->x(1) = x_tet[1];
322 nod_pt->x(2) = x_tet[2];
323 break;
324 case 6:
325 // label 12 in initial sketch: Mid face node on face
326 // spanned by tet nodes 0,2,3
327 s_tet[0] = 1.0 / 3.0;
328 s_tet[1] = 0.0;
329 s_tet[2] = 1.0 / 3.0;
330 nod_pt->set_value(0, s_tet[0]);
331 nod_pt->set_value(1, s_tet[1]);
332 nod_pt->set_value(2, s_tet[2]);
334 nod_pt->x(0) = x_tet[0];
335 nod_pt->x(1) = x_tet[1];
336 nod_pt->x(2) = x_tet[2];
337 break;
338 case 7:
339 // label 14 in initial sketch: Centroid
340 s_tet[0] = 0.25;
341 s_tet[1] = 0.25;
342 s_tet[2] = 0.25;
343 nod_pt->set_value(0, s_tet[0]);
344 nod_pt->set_value(1, s_tet[1]);
345 nod_pt->set_value(2, s_tet[2]);
347 nod_pt->x(0) = x_tet[0];
348 nod_pt->x(1) = x_tet[1];
349 nod_pt->x(2) = x_tet[2];
350 break;
351 }
352 }
353
354 // Create actual zeroth brick element
355 FiniteElement* el_pt = new ELEMENT;
357 Element_pt.push_back(el_pt);
358
359 TFace face0(
361
362 TFace face1(
364
365 TFace face2(
367
368 // Tet vertex nodes along edges emanating from node 0 in brick
370 tet_edge_node[0].resize(2);
371 tet_edge_node[0][0] = 4;
372 tet_edge_node[0][1] = 1;
373 tet_edge_node[1].resize(2);
374 tet_edge_node[1][0] = 6;
375 tet_edge_node[1][1] = 3;
376 tet_edge_node[2].resize(2);
377 tet_edge_node[2][0] = 5;
378 tet_edge_node[2][1] = 2;
379
380 // Node number of tet vertex that node 0 in brick is centred on
381 unsigned central_tet_vertex = 0;
382
383 Node* tet_node_pt = 0;
384 Node* old_node_pt = 0;
385
386 // Corner node
387 {
388 unsigned j = 0;
389
390 // Need new node?
393 if (old_node_pt == 0)
394 {
395 Node* new_node_pt = 0;
396 if (tet_node_pt->is_on_boundary())
397 {
399 }
400 else
401 {
403 }
405 Node_pt.push_back(new_node_pt);
406 Vector<double> s(3);
410 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
412 new_node_pt->x(0) = x_tet[0];
413 new_node_pt->x(1) = x_tet[1];
414 new_node_pt->x(2) = x_tet[2];
415 }
416 // Node already exists
417 else
418 {
420 }
421 }
422
423 // Brick vertex node coindides with mid-edge node on tet edge 0
424 {
425 unsigned j = 2;
426
427 // Need new node?
430 if (old_node_pt == 0)
431 {
432 Node* new_node_pt = 0;
433 if (tet_node_pt->is_on_boundary())
434 {
436 }
437 else
438 {
440 }
442 Node_pt.push_back(new_node_pt);
443 Vector<double> s(3);
447 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
449 new_node_pt->x(0) = x_tet[0];
450 new_node_pt->x(1) = x_tet[1];
451 new_node_pt->x(2) = x_tet[2];
452 }
453 // Node already exists
454 else
455 {
457 }
458 }
459
460 // Brick vertex node coindides with mid vertex node of tet edge 1
461 {
462 unsigned j = 6;
463
464 // Need new node?
467 if (old_node_pt == 0)
468 {
469 Node* new_node_pt = 0;
470 if (tet_node_pt->is_on_boundary())
471 {
473 }
474 else
475 {
477 }
479 Node_pt.push_back(new_node_pt);
480 Vector<double> s(3);
484 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
486 new_node_pt->x(0) = x_tet[0];
487 new_node_pt->x(1) = x_tet[1];
488 new_node_pt->x(2) = x_tet[2];
489 }
490 // Node already exists
491 else
492 {
494 }
495 }
496
497 // Brick vertex node coindides with mid-vertex node of tet edge 2
498 {
499 unsigned j = 18;
500
501 // Need new node?
504 if (old_node_pt == 0)
505 {
506 Node* new_node_pt = 0;
507 if (tet_node_pt->is_on_boundary())
508 {
510 }
511 else
512 {
514 }
516 Node_pt.push_back(new_node_pt);
517 Vector<double> s(3);
521 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
523 new_node_pt->x(0) = x_tet[0];
524 new_node_pt->x(1) = x_tet[1];
525 new_node_pt->x(2) = x_tet[2];
526 }
527 // Node already exists
528 else
529 {
531 }
532 }
533
534 // Brick vertex node in the middle of tet face0, spanned by
535 // tet vertices 0, 1, 2. Enumerated "11" in initial sketch.
536 {
537 unsigned j = 20;
538
539 // Need new node?
541 if (old_node_pt == 0)
542 {
543 Node* new_node_pt = 0;
544 if (face0.is_boundary_face())
545 {
547 }
548 else
549 {
551 }
553 Node_pt.push_back(new_node_pt);
554 Vector<double> s(3);
558 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
560 new_node_pt->x(0) = x_tet[0];
561 new_node_pt->x(1) = x_tet[1];
562 new_node_pt->x(2) = x_tet[2];
563 }
564 // Node already exists
565 else
566 {
568 }
569 }
570
571 // Brick vertex node in the middle of tet face1, spanned by
572 // tet vertices 0, 2, 3. Enumerated "12" in initial sketch.
573 {
574 unsigned j = 24;
575
576 // Need new node?
578 if (old_node_pt == 0)
579 {
580 Node* new_node_pt = 0;
581 if (face1.is_boundary_face())
582 {
584 }
585 else
586 {
588 }
590 Node_pt.push_back(new_node_pt);
591 Vector<double> s(3);
595 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
597 new_node_pt->x(0) = x_tet[0];
598 new_node_pt->x(1) = x_tet[1];
599 new_node_pt->x(2) = x_tet[2];
600 }
601 // Node already exists
602 else
603 {
605 }
606 }
607
608 // Brick vertex node in the middle of tet face2, spanned by
609 // tet vertices 0, 1, 3. Enumerated "13" in initial sketch.
610 {
611 unsigned j = 8;
612
613 // Need new node?
615 if (old_node_pt == 0)
616 {
617 Node* new_node_pt = 0;
618 if (face2.is_boundary_face())
619 {
621 }
622 else
623 {
625 }
627 Node_pt.push_back(new_node_pt);
628 Vector<double> s(3);
632 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
634 new_node_pt->x(0) = x_tet[0];
635 new_node_pt->x(1) = x_tet[1];
636 new_node_pt->x(2) = x_tet[2];
637 }
638 // Node already exists
639 else
640 {
642 }
643 }
644
645 // Brick vertex node in centroid of tet. Only built for first element.
646 // Enumerated "13" in initial sketch.
647 {
648 unsigned j = 26;
649
650 // Always new
651 {
654 Node_pt.push_back(new_node_pt);
655 Vector<double> s(3);
659 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
661 new_node_pt->x(0) = x_tet[0];
662 new_node_pt->x(1) = x_tet[1];
663 new_node_pt->x(2) = x_tet[2];
664 }
665 }
666
667 // Internal brick node -- always built
668 {
669 unsigned j = 13;
670
671 // Always new
672 {
674 Node_pt.push_back(new_node_pt);
675 Vector<double> s(3);
679 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
681 new_node_pt->x(0) = x_tet[0];
682 new_node_pt->x(1) = x_tet[1];
683 new_node_pt->x(2) = x_tet[2];
684 }
685 }
686
687 // Brick edge node between brick nodes 0 and 2
688 {
689 unsigned j = 1;
690
691 // Need new node?
694 if (old_node_pt == 0)
695 {
696 Node* new_node_pt = 0;
697 if (edge.is_boundary_edge())
698 {
700 }
701 else
702 {
704 }
706 Node_pt.push_back(new_node_pt);
707 Vector<double> s(3);
711 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
713 new_node_pt->x(0) = x_tet[0];
714 new_node_pt->x(1) = x_tet[1];
715 new_node_pt->x(2) = x_tet[2];
716 }
717 // Node already exists
718 else
719 {
721 }
722 }
723
724 // Brick edge node between brick nodes 0 and 6
725 {
726 unsigned j = 3;
727
728 // Need new node?
731 if (old_node_pt == 0)
732 {
733 Node* new_node_pt = 0;
734 if (edge.is_boundary_edge())
735 {
737 }
738 else
739 {
741 }
743 Node_pt.push_back(new_node_pt);
744 Vector<double> s(3);
748 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
750 new_node_pt->x(0) = x_tet[0];
751 new_node_pt->x(1) = x_tet[1];
752 new_node_pt->x(2) = x_tet[2];
753 }
754 // Node already exists
755 else
756 {
758 }
759 }
760
761 // Brick edge node between brick nodes 2 and 8
762 {
763 unsigned j = 5;
764
765 // Need new node?
768 if (old_node_pt == 0)
769 {
770 Node* new_node_pt = 0;
771 if (edge.is_boundary_edge())
772 {
774 }
775 else
776 {
778 }
780 Node_pt.push_back(new_node_pt);
781 Vector<double> s(3);
785 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
787 new_node_pt->x(0) = x_tet[0];
788 new_node_pt->x(1) = x_tet[1];
789 new_node_pt->x(2) = x_tet[2];
790 }
791 // Node already exists
792 else
793 {
795 }
796 }
797
798 // Brick edge node between brick nodes 6 and 8
799 {
800 unsigned j = 7;
801
802 // Need new node?
805 if (old_node_pt == 0)
806 {
807 Node* new_node_pt = 0;
808 if (edge.is_boundary_edge())
809 {
811 }
812 else
813 {
815 }
817 Node_pt.push_back(new_node_pt);
818 Vector<double> s(3);
822 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
824 new_node_pt->x(0) = x_tet[0];
825 new_node_pt->x(1) = x_tet[1];
826 new_node_pt->x(2) = x_tet[2];
827 }
828 // Node already exists
829 else
830 {
832 }
833 }
834
835 // Brick edge node between brick nodes 18 and 20
836 {
837 unsigned j = 19;
838
839 // Need new node?
840 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
842 if (old_node_pt == 0)
843 {
844 Node* new_node_pt = 0;
845 if (edge.is_boundary_edge())
846 {
848 }
849 else
850 {
852 }
854 Node_pt.push_back(new_node_pt);
855 Vector<double> s(3);
859 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
861 new_node_pt->x(0) = x_tet[0];
862 new_node_pt->x(1) = x_tet[1];
863 new_node_pt->x(2) = x_tet[2];
864 }
865 // Node already exists
866 else
867 {
869 }
870 }
871
872 // Brick edge node between brick nodes 18 and 24
873 {
874 unsigned j = 21;
875
876 // Need new node?
877 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
879 if (old_node_pt == 0)
880 {
881 Node* new_node_pt = 0;
882 if (edge.is_boundary_edge())
883 {
885 }
886 else
887 {
889 }
891 Node_pt.push_back(new_node_pt);
892 Vector<double> s(3);
896 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
898 new_node_pt->x(0) = x_tet[0];
899 new_node_pt->x(1) = x_tet[1];
900 new_node_pt->x(2) = x_tet[2];
901 }
902 // Node already exists
903 else
904 {
906 }
907 }
908
909 // Brick edge node between brick nodes 20 and 26
910 {
911 unsigned j = 23;
912
913 // Need new node?
914 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
916 if (old_node_pt == 0)
917 {
918 Node* new_node_pt = 0;
919 if (edge.is_boundary_edge())
920 {
922 }
923 else
924 {
926 }
928 Node_pt.push_back(new_node_pt);
929 Vector<double> s(3);
933 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
935 new_node_pt->x(0) = x_tet[0];
936 new_node_pt->x(1) = x_tet[1];
937 new_node_pt->x(2) = x_tet[2];
938 }
939 // Node already exists
940 else
941 {
943 }
944 }
945
946 // Brick edge node between brick nodes 24 and 26
947 {
948 unsigned j = 25;
949
950 // Need new node?
951 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
953 if (old_node_pt == 0)
954 {
955 Node* new_node_pt = 0;
956 if (edge.is_boundary_edge())
957 {
959 }
960 else
961 {
963 }
965 Node_pt.push_back(new_node_pt);
966 Vector<double> s(3);
970 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
972 new_node_pt->x(0) = x_tet[0];
973 new_node_pt->x(1) = x_tet[1];
974 new_node_pt->x(2) = x_tet[2];
975 }
976 // Node already exists
977 else
978 {
980 }
981 }
982
983 // Brick edge node between brick nodes 0 and 18
984 {
985 unsigned j = 9;
986
987 // Need new node?
988 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
990 if (old_node_pt == 0)
991 {
992 Node* new_node_pt = 0;
993 if (edge.is_boundary_edge())
994 {
996 }
997 else
998 {
1000 }
1002 Node_pt.push_back(new_node_pt);
1003 Vector<double> s(3);
1007 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1009 new_node_pt->x(0) = x_tet[0];
1010 new_node_pt->x(1) = x_tet[1];
1011 new_node_pt->x(2) = x_tet[2];
1012 }
1013 // Node already exists
1014 else
1015 {
1017 }
1018 }
1019
1020 // Brick edge node between brick nodes 2 and 20
1021 {
1022 unsigned j = 11;
1023
1024 // Need new node?
1025 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
1027 if (old_node_pt == 0)
1028 {
1029 Node* new_node_pt = 0;
1030 if (edge.is_boundary_edge())
1031 {
1033 }
1034 else
1035 {
1037 }
1039 Node_pt.push_back(new_node_pt);
1040 Vector<double> s(3);
1044 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1046 new_node_pt->x(0) = x_tet[0];
1047 new_node_pt->x(1) = x_tet[1];
1048 new_node_pt->x(2) = x_tet[2];
1049 }
1050 // Node already exists
1051 else
1052 {
1054 }
1055 }
1056
1057 // Brick edge node between brick nodes 6 and 24
1058 {
1059 unsigned j = 15;
1060
1061 // Need new node?
1062 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
1064 if (old_node_pt == 0)
1065 {
1066 Node* new_node_pt = 0;
1067 if (edge.is_boundary_edge())
1068 {
1070 }
1071 else
1072 {
1074 }
1076 Node_pt.push_back(new_node_pt);
1077 Vector<double> s(3);
1081 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1083 new_node_pt->x(0) = x_tet[0];
1084 new_node_pt->x(1) = x_tet[1];
1085 new_node_pt->x(2) = x_tet[2];
1086 }
1087 // Node already exists
1088 else
1089 {
1091 }
1092 }
1093
1094 // Brick edge node between brick nodes 8 and 26
1095 {
1096 unsigned j = 17;
1097
1098 // Need new node?
1099 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
1101 if (old_node_pt == 0)
1102 {
1103 Node* new_node_pt = 0;
1104 if (edge.is_boundary_edge())
1105 {
1107 }
1108 else
1109 {
1111 }
1113 Node_pt.push_back(new_node_pt);
1114 Vector<double> s(3);
1118 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1120 new_node_pt->x(0) = x_tet[0];
1121 new_node_pt->x(1) = x_tet[1];
1122 new_node_pt->x(2) = x_tet[2];
1123 }
1124 // Node already exists
1125 else
1126 {
1128 }
1129 }
1130
1131 // Mid brick-face node associated with face
1132 // spanned by mid-vertex nodes associated with tet edges 0 and 2
1133 {
1134 unsigned j = 10;
1135
1136 // Need new node?
1140
1142 if (old_node_pt == 0)
1143 {
1144 Node* new_node_pt = 0;
1145 if (face.is_boundary_face())
1146 {
1148 }
1149 else
1150 {
1152 }
1154 Node_pt.push_back(new_node_pt);
1155 Vector<double> s(3);
1159 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1161 new_node_pt->x(0) = x_tet[0];
1162 new_node_pt->x(1) = x_tet[1];
1163 new_node_pt->x(2) = x_tet[2];
1164 }
1165 // Node already exists
1166 else
1167 {
1169 }
1170 }
1171
1172 // Mid brick-face node associated with face
1173 // spanned by mid-vertex nodes associated with tet edges 1 and 2
1174 {
1175 unsigned j = 12;
1176
1177 // Need new node?
1181
1183 if (old_node_pt == 0)
1184 {
1185 Node* new_node_pt = 0;
1186 if (face.is_boundary_face())
1187 {
1189 }
1190 else
1191 {
1193 }
1195 Node_pt.push_back(new_node_pt);
1196 Vector<double> s(3);
1200 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1202 new_node_pt->x(0) = x_tet[0];
1203 new_node_pt->x(1) = x_tet[1];
1204 new_node_pt->x(2) = x_tet[2];
1205 }
1206 // Node already exists
1207 else
1208 {
1210 }
1211 }
1212
1213 // Mid brick-face node associated with face
1214 // spanned by mid-vertex nodes associated with tet edges 0 and 1
1215 {
1216 unsigned j = 4;
1217
1218 // Need new node?
1222
1224 if (old_node_pt == 0)
1225 {
1226 Node* new_node_pt = 0;
1227 if (face.is_boundary_face())
1228 {
1230 }
1231 else
1232 {
1234 }
1236 Node_pt.push_back(new_node_pt);
1237 Vector<double> s(3);
1241 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1243 new_node_pt->x(0) = x_tet[0];
1244 new_node_pt->x(1) = x_tet[1];
1245 new_node_pt->x(2) = x_tet[2];
1246 }
1247 // Node already exists
1248 else
1249 {
1251 }
1252 }
1253
1254 // Top mid brick-face node -- only built by first element
1255 {
1256 unsigned j = 22;
1258 Node_pt.push_back(new_node_pt);
1259 Vector<double> s(3);
1263 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1266 new_node_pt->x(0) = x_tet[0];
1267 new_node_pt->x(1) = x_tet[1];
1268 new_node_pt->x(2) = x_tet[2];
1269 }
1270
1271 // Right mid brick-face node -- only built by first element
1272 {
1273 unsigned j = 14;
1275 Node_pt.push_back(new_node_pt);
1276 Vector<double> s(3);
1280 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1283 new_node_pt->x(0) = x_tet[0];
1284 new_node_pt->x(1) = x_tet[1];
1285 new_node_pt->x(2) = x_tet[2];
1286 }
1287
1288 // Back mid brick-face node -- only built by first element
1289 {
1290 unsigned j = 16;
1292 Node_pt.push_back(new_node_pt);
1293 Vector<double> s(3);
1297 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
1300 new_node_pt->x(0) = x_tet[0];
1301 new_node_pt->x(1) = x_tet[1];
1302 new_node_pt->x(2) = x_tet[2];
1303 }
1304 }
1305
1306 // Second brick element is centred at node 1 of tet:
1307 //--------------------------------------------------
1308 {
1309 // Assign coordinates of dummy element
1310 for (unsigned j = 0; j < 8; j++)
1311 {
1315 switch (j)
1316 {
1317 case 0:
1319 nod_pt->set_value(0, s_tet[0]);
1320 nod_pt->set_value(1, s_tet[1]);
1321 nod_pt->set_value(2, s_tet[2]);
1323 nod_pt->x(0) = x_tet[0];
1324 nod_pt->x(1) = x_tet[1];
1325 nod_pt->x(2) = x_tet[2];
1326 break;
1327 case 1:
1329 nod_pt->set_value(0, s_tet[0]);
1330 nod_pt->set_value(1, s_tet[1]);
1331 nod_pt->set_value(2, s_tet[2]);
1333 nod_pt->x(0) = x_tet[0];
1334 nod_pt->x(1) = x_tet[1];
1335 nod_pt->x(2) = x_tet[2];
1336 break;
1337 case 2:
1339 nod_pt->set_value(0, s_tet[0]);
1340 nod_pt->set_value(1, s_tet[1]);
1341 nod_pt->set_value(2, s_tet[2]);
1343 nod_pt->x(0) = x_tet[0];
1344 nod_pt->x(1) = x_tet[1];
1345 nod_pt->x(2) = x_tet[2];
1346 break;
1347 case 3:
1348 // label 13 in initial sketch: Mid face node on face
1349 // spanned by tet nodes 0,1,3
1350 s_tet[0] = 1.0 / 3.0;
1351 s_tet[1] = 1.0 / 3.0;
1352 s_tet[2] = 0.0;
1353 nod_pt->set_value(0, s_tet[0]);
1354 nod_pt->set_value(1, s_tet[1]);
1355 nod_pt->set_value(2, s_tet[2]);
1357 nod_pt->x(0) = x_tet[0];
1358 nod_pt->x(1) = x_tet[1];
1359 nod_pt->x(2) = x_tet[2];
1360 break;
1361 case 4:
1363 nod_pt->set_value(0, s_tet[0]);
1364 nod_pt->set_value(1, s_tet[1]);
1365 nod_pt->set_value(2, s_tet[2]);
1367 nod_pt->x(0) = x_tet[0];
1368 nod_pt->x(1) = x_tet[1];
1369 nod_pt->x(2) = x_tet[2];
1370 break;
1371 case 5:
1372 // label 10 in initial sketch: Mid face node on face
1373 // spanned by tet nodes 1,2,3
1374 s_tet[0] = 0.0;
1375 s_tet[1] = 1.0 / 3.0;
1376 s_tet[2] = 1.0 / 3.0;
1377 nod_pt->set_value(0, s_tet[0]);
1378 nod_pt->set_value(1, s_tet[1]);
1379 nod_pt->set_value(2, s_tet[2]);
1381 nod_pt->x(0) = x_tet[0];
1382 nod_pt->x(1) = x_tet[1];
1383 nod_pt->x(2) = x_tet[2];
1384 break;
1385 case 6:
1386 // label 11 in initial sketch: Mid face node on face
1387 // spanned by tet nodes 0,1,2
1388 s_tet[0] = 1.0 / 3.0;
1389 s_tet[1] = 1.0 / 3.0;
1390 s_tet[2] = 1.0 / 3.0;
1391 nod_pt->set_value(0, s_tet[0]);
1392 nod_pt->set_value(1, s_tet[1]);
1393 nod_pt->set_value(2, s_tet[2]);
1395 nod_pt->x(0) = x_tet[0];
1396 nod_pt->x(1) = x_tet[1];
1397 nod_pt->x(2) = x_tet[2];
1398 break;
1399 case 7:
1400 // label 14 in initial sketch: Centroid
1401 s_tet[0] = 0.25;
1402 s_tet[1] = 0.25;
1403 s_tet[2] = 0.25;
1404 nod_pt->set_value(0, s_tet[0]);
1405 nod_pt->set_value(1, s_tet[1]);
1406 nod_pt->set_value(2, s_tet[2]);
1408 nod_pt->x(0) = x_tet[0];
1409 nod_pt->x(1) = x_tet[1];
1410 nod_pt->x(2) = x_tet[2];
1411 break;
1412 }
1413 }
1414
1415 // Create actual first brick element
1416 FiniteElement* el_pt = new ELEMENT;
1418 Element_pt.push_back(el_pt);
1419
1420 TFace face0(
1422
1423 TFace face1(
1425
1426 TFace face2(
1428
1429 // Tet vertex nodes along edges emanating from node 0 in brick
1431 tet_edge_node[0].resize(2);
1432 tet_edge_node[0][0] = 9;
1433 tet_edge_node[0][1] = 3;
1434 tet_edge_node[1].resize(2);
1435 tet_edge_node[1][0] = 4;
1436 tet_edge_node[1][1] = 0;
1437 tet_edge_node[2].resize(2);
1438 tet_edge_node[2][0] = 7;
1439 tet_edge_node[2][1] = 2;
1440
1441 // Node number of tet vertex that node 0 in brick is centred on
1442 unsigned central_tet_vertex = 1;
1443
1444 Node* tet_node_pt = 0;
1445 Node* old_node_pt = 0;
1446
1447 // Corner node
1448 {
1449 unsigned j = 0;
1450
1451 // Need new node?
1454 if (old_node_pt == 0)
1455 {
1456 Node* new_node_pt = 0;
1457 if (tet_node_pt->is_on_boundary())
1458 {
1460 }
1461 else
1462 {
1464 }
1466 Node_pt.push_back(new_node_pt);
1467 Vector<double> s(3);
1471 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1473 new_node_pt->x(0) = x_tet[0];
1474 new_node_pt->x(1) = x_tet[1];
1475 new_node_pt->x(2) = x_tet[2];
1476 }
1477 // Node already exists
1478 else
1479 {
1481 }
1482 }
1483
1484 // Brick vertex node coindides with mid-edge node on tet edge 0
1485 {
1486 unsigned j = 2;
1487
1488 // Need new node?
1491 if (old_node_pt == 0)
1492 {
1493 Node* new_node_pt = 0;
1494 if (tet_node_pt->is_on_boundary())
1495 {
1497 }
1498 else
1499 {
1501 }
1503 Node_pt.push_back(new_node_pt);
1504 Vector<double> s(3);
1508 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1510 new_node_pt->x(0) = x_tet[0];
1511 new_node_pt->x(1) = x_tet[1];
1512 new_node_pt->x(2) = x_tet[2];
1513 }
1514 // Node already exists
1515 else
1516 {
1518 }
1519 }
1520
1521 // Brick vertex node coindides with mid vertex node of tet edge 1
1522 {
1523 unsigned j = 6;
1524
1525 // Need new node?
1528 if (old_node_pt == 0)
1529 {
1530 Node* new_node_pt = 0;
1531 if (tet_node_pt->is_on_boundary())
1532 {
1534 }
1535 else
1536 {
1538 }
1540 Node_pt.push_back(new_node_pt);
1541 Vector<double> s(3);
1545 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1547 new_node_pt->x(0) = x_tet[0];
1548 new_node_pt->x(1) = x_tet[1];
1549 new_node_pt->x(2) = x_tet[2];
1550 }
1551 // Node already exists
1552 else
1553 {
1555 }
1556 }
1557
1558 // Brick vertex node coindides with mid-vertex node of tet edge 2
1559 {
1560 unsigned j = 18;
1561
1562 // Need new node?
1565 if (old_node_pt == 0)
1566 {
1567 Node* new_node_pt = 0;
1568 if (tet_node_pt->is_on_boundary())
1569 {
1571 }
1572 else
1573 {
1575 }
1577 Node_pt.push_back(new_node_pt);
1578 Vector<double> s(3);
1582 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1584 new_node_pt->x(0) = x_tet[0];
1585 new_node_pt->x(1) = x_tet[1];
1586 new_node_pt->x(2) = x_tet[2];
1587 }
1588 // Node already exists
1589 else
1590 {
1592 }
1593 }
1594
1595 // Brick vertex node in the middle of tet face0
1596 {
1597 unsigned j = 20;
1598
1599 // Need new node?
1601 if (old_node_pt == 0)
1602 {
1603 Node* new_node_pt = 0;
1604 if (face0.is_boundary_face())
1605 {
1607 }
1608 else
1609 {
1611 }
1613 Node_pt.push_back(new_node_pt);
1614 Vector<double> s(3);
1618 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1620 new_node_pt->x(0) = x_tet[0];
1621 new_node_pt->x(1) = x_tet[1];
1622 new_node_pt->x(2) = x_tet[2];
1623 }
1624 // Node already exists
1625 else
1626 {
1628 }
1629 }
1630
1631 // Brick vertex node in the middle of tet face1
1632 {
1633 unsigned j = 24;
1634
1635 // Need new node?
1637 if (old_node_pt == 0)
1638 {
1639 Node* new_node_pt = 0;
1640 if (face1.is_boundary_face())
1641 {
1643 }
1644 else
1645 {
1647 }
1649 Node_pt.push_back(new_node_pt);
1650 Vector<double> s(3);
1654 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1656 new_node_pt->x(0) = x_tet[0];
1657 new_node_pt->x(1) = x_tet[1];
1658 new_node_pt->x(2) = x_tet[2];
1659 }
1660 // Node already exists
1661 else
1662 {
1664 }
1665 }
1666
1667 // Brick vertex node in the middle of face2
1668 {
1669 unsigned j = 8;
1670
1671 // Need new node?
1673 if (old_node_pt == 0)
1674 {
1675 Node* new_node_pt = 0;
1676 if (face2.is_boundary_face())
1677 {
1679 }
1680 else
1681 {
1683 }
1685 Node_pt.push_back(new_node_pt);
1686 Vector<double> s(3);
1690 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1692 new_node_pt->x(0) = x_tet[0];
1693 new_node_pt->x(1) = x_tet[1];
1694 new_node_pt->x(2) = x_tet[2];
1695 }
1696 // Node already exists
1697 else
1698 {
1700 }
1701 }
1702
1703 // Brick vertex node in centroid of tet. Only built for first element.
1704 // Enumerated "13" in initial sketch.
1705 {
1706 unsigned j = 26;
1707
1708 // Always copied
1710 }
1711
1712 // Internal brick node -- always built
1713 {
1714 unsigned j = 13;
1715
1716 // Always new
1717 {
1719 Node_pt.push_back(new_node_pt);
1720 Vector<double> s(3);
1724 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1726 new_node_pt->x(0) = x_tet[0];
1727 new_node_pt->x(1) = x_tet[1];
1728 new_node_pt->x(2) = x_tet[2];
1729 }
1730 }
1731
1732 // Brick edge node between brick nodes 0 and 2
1733 {
1734 unsigned j = 1;
1735
1736 // Need new node?
1737 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
1739 if (old_node_pt == 0)
1740 {
1741 Node* new_node_pt = 0;
1742 if (edge.is_boundary_edge())
1743 {
1745 }
1746 else
1747 {
1749 }
1751 Node_pt.push_back(new_node_pt);
1752 Vector<double> s(3);
1756 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1758 new_node_pt->x(0) = x_tet[0];
1759 new_node_pt->x(1) = x_tet[1];
1760 new_node_pt->x(2) = x_tet[2];
1761 }
1762 // Node already exists
1763 else
1764 {
1766 }
1767 }
1768
1769 // Brick edge node between brick nodes 0 and 6
1770 {
1771 unsigned j = 3;
1772
1773 // Need new node?
1774 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
1776 if (old_node_pt == 0)
1777 {
1778 Node* new_node_pt = 0;
1779 if (edge.is_boundary_edge())
1780 {
1782 }
1783 else
1784 {
1786 }
1788 Node_pt.push_back(new_node_pt);
1789 Vector<double> s(3);
1793 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1795 new_node_pt->x(0) = x_tet[0];
1796 new_node_pt->x(1) = x_tet[1];
1797 new_node_pt->x(2) = x_tet[2];
1798 }
1799 // Node already exists
1800 else
1801 {
1803 }
1804 }
1805
1806 // Brick edge node between brick nodes 2 and 8
1807 {
1808 unsigned j = 5;
1809
1810 // Need new node?
1811 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
1813 if (old_node_pt == 0)
1814 {
1815 Node* new_node_pt = 0;
1816 if (edge.is_boundary_edge())
1817 {
1819 }
1820 else
1821 {
1823 }
1825 Node_pt.push_back(new_node_pt);
1826 Vector<double> s(3);
1830 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1832 new_node_pt->x(0) = x_tet[0];
1833 new_node_pt->x(1) = x_tet[1];
1834 new_node_pt->x(2) = x_tet[2];
1835 }
1836 // Node already exists
1837 else
1838 {
1840 }
1841 }
1842
1843 // Brick edge node between brick nodes 6 and 8
1844 {
1845 unsigned j = 7;
1846
1847 // Need new node?
1848 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
1850 if (old_node_pt == 0)
1851 {
1852 Node* new_node_pt = 0;
1853 if (edge.is_boundary_edge())
1854 {
1856 }
1857 else
1858 {
1860 }
1862 Node_pt.push_back(new_node_pt);
1863 Vector<double> s(3);
1867 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1869 new_node_pt->x(0) = x_tet[0];
1870 new_node_pt->x(1) = x_tet[1];
1871 new_node_pt->x(2) = x_tet[2];
1872 }
1873 // Node already exists
1874 else
1875 {
1877 }
1878 }
1879
1880 // Brick edge node between brick nodes 18 and 20
1881 {
1882 unsigned j = 19;
1883
1884 // Need new node?
1885 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
1887 if (old_node_pt == 0)
1888 {
1889 Node* new_node_pt = 0;
1890 if (edge.is_boundary_edge())
1891 {
1893 }
1894 else
1895 {
1897 }
1899 Node_pt.push_back(new_node_pt);
1900 Vector<double> s(3);
1904 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1906 new_node_pt->x(0) = x_tet[0];
1907 new_node_pt->x(1) = x_tet[1];
1908 new_node_pt->x(2) = x_tet[2];
1909 }
1910 // Node already exists
1911 else
1912 {
1914 }
1915 }
1916
1917 // Brick edge node between brick nodes 18 and 24
1918 {
1919 unsigned j = 21;
1920
1921 // Need new node?
1922 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
1924 if (old_node_pt == 0)
1925 {
1926 Node* new_node_pt = 0;
1927 if (edge.is_boundary_edge())
1928 {
1930 }
1931 else
1932 {
1934 }
1936 Node_pt.push_back(new_node_pt);
1937 Vector<double> s(3);
1941 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1943 new_node_pt->x(0) = x_tet[0];
1944 new_node_pt->x(1) = x_tet[1];
1945 new_node_pt->x(2) = x_tet[2];
1946 }
1947 // Node already exists
1948 else
1949 {
1951 }
1952 }
1953
1954 // Brick edge node between brick nodes 20 and 26
1955 {
1956 unsigned j = 23;
1957
1958 // Need new node?
1959 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
1961 if (old_node_pt == 0)
1962 {
1963 Node* new_node_pt = 0;
1964 if (edge.is_boundary_edge())
1965 {
1967 }
1968 else
1969 {
1971 }
1973 Node_pt.push_back(new_node_pt);
1974 Vector<double> s(3);
1978 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
1980 new_node_pt->x(0) = x_tet[0];
1981 new_node_pt->x(1) = x_tet[1];
1982 new_node_pt->x(2) = x_tet[2];
1983 }
1984 // Node already exists
1985 else
1986 {
1988 }
1989 }
1990
1991 // Brick edge node between brick nodes 24 and 26
1992 {
1993 unsigned j = 25;
1994
1995 // Need new node?
1996 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
1998 if (old_node_pt == 0)
1999 {
2000 Node* new_node_pt = 0;
2001 if (edge.is_boundary_edge())
2002 {
2004 }
2005 else
2006 {
2008 }
2010 Node_pt.push_back(new_node_pt);
2011 Vector<double> s(3);
2015 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2017 new_node_pt->x(0) = x_tet[0];
2018 new_node_pt->x(1) = x_tet[1];
2019 new_node_pt->x(2) = x_tet[2];
2020 }
2021 // Node already exists
2022 else
2023 {
2025 }
2026 }
2027
2028 // Brick edge node between brick nodes 0 and 18
2029 {
2030 unsigned j = 9;
2031
2032 // Need new node?
2033 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
2035 if (old_node_pt == 0)
2036 {
2037 Node* new_node_pt = 0;
2038 if (edge.is_boundary_edge())
2039 {
2041 }
2042 else
2043 {
2045 }
2047 Node_pt.push_back(new_node_pt);
2048 Vector<double> s(3);
2052 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2054 new_node_pt->x(0) = x_tet[0];
2055 new_node_pt->x(1) = x_tet[1];
2056 new_node_pt->x(2) = x_tet[2];
2057 }
2058 // Node already exists
2059 else
2060 {
2062 }
2063 }
2064
2065 // Brick edge node between brick nodes 2 and 20
2066 {
2067 unsigned j = 11;
2068
2069 // Need new node?
2070 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
2072 if (old_node_pt == 0)
2073 {
2074 Node* new_node_pt = 0;
2075 if (edge.is_boundary_edge())
2076 {
2078 }
2079 else
2080 {
2082 }
2084 Node_pt.push_back(new_node_pt);
2085 Vector<double> s(3);
2089 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2091 new_node_pt->x(0) = x_tet[0];
2092 new_node_pt->x(1) = x_tet[1];
2093 new_node_pt->x(2) = x_tet[2];
2094 }
2095 // Node already exists
2096 else
2097 {
2099 }
2100 }
2101
2102 // Brick edge node between brick nodes 6 and 24
2103 {
2104 unsigned j = 15;
2105
2106 // Need new node?
2107 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
2109 if (old_node_pt == 0)
2110 {
2111 Node* new_node_pt = 0;
2112 if (edge.is_boundary_edge())
2113 {
2115 }
2116 else
2117 {
2119 }
2121 Node_pt.push_back(new_node_pt);
2122 Vector<double> s(3);
2126 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2128 new_node_pt->x(0) = x_tet[0];
2129 new_node_pt->x(1) = x_tet[1];
2130 new_node_pt->x(2) = x_tet[2];
2131 }
2132 // Node already exists
2133 else
2134 {
2136 }
2137 }
2138
2139 // Brick edge node between brick nodes 8 and 26
2140 {
2141 unsigned j = 17;
2142
2143 // Need new node?
2144 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
2146 if (old_node_pt == 0)
2147 {
2148 Node* new_node_pt = 0;
2149 if (edge.is_boundary_edge())
2150 {
2152 }
2153 else
2154 {
2156 }
2158 Node_pt.push_back(new_node_pt);
2159 Vector<double> s(3);
2163 dummy_q_el_pt[1]->interpolated_s_tet(s, s_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 {
2173 }
2174 }
2175
2176 // Mid brick-face node associated with face
2177 // spanned by mid-vertex nodes associated with tet edges 0 and 2
2178 {
2179 unsigned j = 10;
2180
2181 // Need new node?
2185
2187 if (old_node_pt == 0)
2188 {
2189 Node* new_node_pt = 0;
2190 if (face.is_boundary_face())
2191 {
2193 }
2194 else
2195 {
2197 }
2199 Node_pt.push_back(new_node_pt);
2200 Vector<double> s(3);
2204 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2206 new_node_pt->x(0) = x_tet[0];
2207 new_node_pt->x(1) = x_tet[1];
2208 new_node_pt->x(2) = x_tet[2];
2209 }
2210 // Node already exists
2211 else
2212 {
2214 }
2215 }
2216
2217 // Mid brick-face node associated with face
2218 // spanned by mid-vertex nodes associated with tet edges 1 and 2
2219 {
2220 unsigned j = 12;
2221
2222 // Need new node?
2226
2228 if (old_node_pt == 0)
2229 {
2230 Node* new_node_pt = 0;
2231 if (face.is_boundary_face())
2232 {
2234 }
2235 else
2236 {
2238 }
2240 Node_pt.push_back(new_node_pt);
2241 Vector<double> s(3);
2245 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2247 new_node_pt->x(0) = x_tet[0];
2248 new_node_pt->x(1) = x_tet[1];
2249 new_node_pt->x(2) = x_tet[2];
2250 }
2251 // Node already exists
2252 else
2253 {
2255 }
2256 }
2257
2258 // Mid brick-face node associated with face
2259 // spanned by mid-vertex nodes associated with tet edges 0 and 1
2260 {
2261 unsigned j = 4;
2262
2263 // Need new node?
2267
2269 if (old_node_pt == 0)
2270 {
2271 Node* new_node_pt = 0;
2272 if (face.is_boundary_face())
2273 {
2275 }
2276 else
2277 {
2279 }
2281 Node_pt.push_back(new_node_pt);
2282 Vector<double> s(3);
2286 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2288 new_node_pt->x(0) = x_tet[0];
2289 new_node_pt->x(1) = x_tet[1];
2290 new_node_pt->x(2) = x_tet[2];
2291 }
2292 // Node already exists
2293 else
2294 {
2296 }
2297 }
2298
2299 // Top mid brick-face node -- only built by this element
2300 {
2301 unsigned j = 22;
2303 Node_pt.push_back(new_node_pt);
2304 Vector<double> s(3);
2308 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2311 new_node_pt->x(0) = x_tet[0];
2312 new_node_pt->x(1) = x_tet[1];
2313 new_node_pt->x(2) = x_tet[2];
2314 }
2315
2316 // Right mid brick-face node -- only built by this element
2317 {
2318 unsigned j = 14;
2320 Node_pt.push_back(new_node_pt);
2321 Vector<double> s(3);
2325 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
2328 new_node_pt->x(0) = x_tet[0];
2329 new_node_pt->x(1) = x_tet[1];
2330 new_node_pt->x(2) = x_tet[2];
2331 }
2332
2333 // Back mid brick-face node copy from previous element
2334 {
2335 unsigned j = 16;
2336
2337 // Always copied
2339 }
2340 }
2341
2342 // Third brick element is centred at node 3 of tet:
2343 //-------------------------------------------------
2344 {
2345 // Assign coordinates of dummy element
2346 for (unsigned j = 0; j < 8; j++)
2347 {
2351 switch (j)
2352 {
2353 case 0:
2355 nod_pt->set_value(0, s_tet[0]);
2356 nod_pt->set_value(1, s_tet[1]);
2357 nod_pt->set_value(2, s_tet[2]);
2359 nod_pt->x(0) = x_tet[0];
2360 nod_pt->x(1) = x_tet[1];
2361 nod_pt->x(2) = x_tet[2];
2362 break;
2363 case 1:
2365 nod_pt->set_value(0, s_tet[0]);
2366 nod_pt->set_value(1, s_tet[1]);
2367 nod_pt->set_value(2, s_tet[2]);
2369 nod_pt->x(0) = x_tet[0];
2370 nod_pt->x(1) = x_tet[1];
2371 nod_pt->x(2) = x_tet[2];
2372 break;
2373 case 2:
2375 nod_pt->set_value(0, s_tet[0]);
2376 nod_pt->set_value(1, s_tet[1]);
2377 nod_pt->set_value(2, s_tet[2]);
2379 nod_pt->x(0) = x_tet[0];
2380 nod_pt->x(1) = x_tet[1];
2381 nod_pt->x(2) = x_tet[2];
2382 break;
2383 case 3:
2384 // label 13 in initial sketch: Mid face node on face
2385 // spanned by tet nodes 0,1,3
2386 s_tet[0] = 1.0 / 3.0;
2387 s_tet[1] = 1.0 / 3.0;
2388 s_tet[2] = 0.0;
2389 nod_pt->set_value(0, s_tet[0]);
2390 nod_pt->set_value(1, s_tet[1]);
2391 nod_pt->set_value(2, s_tet[2]);
2393 nod_pt->x(0) = x_tet[0];
2394 nod_pt->x(1) = x_tet[1];
2395 nod_pt->x(2) = x_tet[2];
2396 break;
2397 case 4:
2399 nod_pt->set_value(0, s_tet[0]);
2400 nod_pt->set_value(1, s_tet[1]);
2401 nod_pt->set_value(2, s_tet[2]);
2403 nod_pt->x(0) = x_tet[0];
2404 nod_pt->x(1) = x_tet[1];
2405 nod_pt->x(2) = x_tet[2];
2406 break;
2407 case 5:
2408 // label 12 in initial sketch: Mid face node on face
2409 // spanned by tet nodes 0,2,3
2410 s_tet[0] = 1.0 / 3.0;
2411 s_tet[1] = 0.0;
2412 s_tet[2] = 1.0 / 3.0;
2413 nod_pt->set_value(0, s_tet[0]);
2414 nod_pt->set_value(1, s_tet[1]);
2415 nod_pt->set_value(2, s_tet[2]);
2417 nod_pt->x(0) = x_tet[0];
2418 nod_pt->x(1) = x_tet[1];
2419 nod_pt->x(2) = x_tet[2];
2420 break;
2421 case 6:
2422 // label 10 in initial sketch: Mid face node on face
2423 // spanned by tet nodes 1,2,3
2424 s_tet[0] = 0.0;
2425 s_tet[1] = 1.0 / 3.0;
2426 s_tet[2] = 1.0 / 3.0;
2427 nod_pt->set_value(0, s_tet[0]);
2428 nod_pt->set_value(1, s_tet[1]);
2429 nod_pt->set_value(2, s_tet[2]);
2431 nod_pt->x(0) = x_tet[0];
2432 nod_pt->x(1) = x_tet[1];
2433 nod_pt->x(2) = x_tet[2];
2434 break;
2435 case 7:
2436 // label 14 in initial sketch: Centroid
2437 s_tet[0] = 0.25;
2438 s_tet[1] = 0.25;
2439 s_tet[2] = 0.25;
2440 nod_pt->set_value(0, s_tet[0]);
2441 nod_pt->set_value(1, s_tet[1]);
2442 nod_pt->set_value(2, s_tet[2]);
2444 nod_pt->x(0) = x_tet[0];
2445 nod_pt->x(1) = x_tet[1];
2446 nod_pt->x(2) = x_tet[2];
2447 break;
2448 }
2449 }
2450
2451 // Create actual second brick element
2452 FiniteElement* el_pt = new ELEMENT;
2454 Element_pt.push_back(el_pt);
2455
2456 TFace face0(
2458
2459 TFace face1(
2461
2462 TFace face2(
2464
2465 // Tet vertex nodes along edges emanating from node 0 in brick
2467 tet_edge_node[0].resize(2);
2468 tet_edge_node[0][0] = 6;
2469 tet_edge_node[0][1] = 0;
2470 tet_edge_node[1].resize(2);
2471 tet_edge_node[1][0] = 9;
2472 tet_edge_node[1][1] = 1;
2473 tet_edge_node[2].resize(2);
2474 tet_edge_node[2][0] = 8;
2475 tet_edge_node[2][1] = 2;
2476
2477 // Node number of tet vertex that node 0 in brick is centred on
2478 unsigned central_tet_vertex = 3;
2479
2480 Node* tet_node_pt = 0;
2481 Node* old_node_pt = 0;
2482
2483 // Corner node
2484 {
2485 unsigned j = 0;
2486
2487 // Need new node?
2490 if (old_node_pt == 0)
2491 {
2492 Node* new_node_pt = 0;
2493 if (tet_node_pt->is_on_boundary())
2494 {
2496 }
2497 else
2498 {
2500 }
2502 Node_pt.push_back(new_node_pt);
2503 Vector<double> s(3);
2507 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2509 new_node_pt->x(0) = x_tet[0];
2510 new_node_pt->x(1) = x_tet[1];
2511 new_node_pt->x(2) = x_tet[2];
2512 }
2513 // Node already exists
2514 else
2515 {
2517 }
2518 }
2519
2520 // Brick vertex node coindides with mid-edge node on tet edge 0
2521 {
2522 unsigned j = 2;
2523
2524 // Need new node?
2527 if (old_node_pt == 0)
2528 {
2529 Node* new_node_pt = 0;
2530 if (tet_node_pt->is_on_boundary())
2531 {
2533 }
2534 else
2535 {
2537 }
2539 Node_pt.push_back(new_node_pt);
2540 Vector<double> s(3);
2544 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2546 new_node_pt->x(0) = x_tet[0];
2547 new_node_pt->x(1) = x_tet[1];
2548 new_node_pt->x(2) = x_tet[2];
2549 }
2550 // Node already exists
2551 else
2552 {
2554 }
2555 }
2556
2557 // Brick vertex node coindides with mid vertex node of tet edge 1
2558 {
2559 unsigned j = 6;
2560
2561 // Need new node?
2564 if (old_node_pt == 0)
2565 {
2566 Node* new_node_pt = 0;
2567 if (tet_node_pt->is_on_boundary())
2568 {
2570 }
2571 else
2572 {
2574 }
2576 Node_pt.push_back(new_node_pt);
2577 Vector<double> s(3);
2581 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2583 new_node_pt->x(0) = x_tet[0];
2584 new_node_pt->x(1) = x_tet[1];
2585 new_node_pt->x(2) = x_tet[2];
2586 }
2587 // Node already exists
2588 else
2589 {
2591 }
2592 }
2593
2594 // Brick vertex node coindides with mid-vertex node of tet edge 2
2595 {
2596 unsigned j = 18;
2597
2598 // Need new node?
2601 if (old_node_pt == 0)
2602 {
2603 Node* new_node_pt = 0;
2604 if (tet_node_pt->is_on_boundary())
2605 {
2607 }
2608 else
2609 {
2611 }
2613 Node_pt.push_back(new_node_pt);
2614 Vector<double> s(3);
2618 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2620 new_node_pt->x(0) = x_tet[0];
2621 new_node_pt->x(1) = x_tet[1];
2622 new_node_pt->x(2) = x_tet[2];
2623 }
2624 // Node already exists
2625 else
2626 {
2628 }
2629 }
2630
2631 // Brick vertex node in the middle of tet face0
2632 {
2633 unsigned j = 20;
2634
2635 // Need new node?
2637 if (old_node_pt == 0)
2638 {
2639 Node* new_node_pt = 0;
2640 if (face0.is_boundary_face())
2641 {
2643 }
2644 else
2645 {
2647 }
2649 Node_pt.push_back(new_node_pt);
2650 Vector<double> s(3);
2654 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2656 new_node_pt->x(0) = x_tet[0];
2657 new_node_pt->x(1) = x_tet[1];
2658 new_node_pt->x(2) = x_tet[2];
2659 }
2660 // Node already exists
2661 else
2662 {
2664 }
2665 }
2666
2667 // Brick vertex node in the middle of tet face1
2668 {
2669 unsigned j = 24;
2670
2671 // Need new node?
2673 if (old_node_pt == 0)
2674 {
2675 Node* new_node_pt = 0;
2676 if (face1.is_boundary_face())
2677 {
2679 }
2680 else
2681 {
2683 }
2685 Node_pt.push_back(new_node_pt);
2686 Vector<double> s(3);
2690 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2692 new_node_pt->x(0) = x_tet[0];
2693 new_node_pt->x(1) = x_tet[1];
2694 new_node_pt->x(2) = x_tet[2];
2695 }
2696 // Node already exists
2697 else
2698 {
2700 }
2701 }
2702
2703 // Brick vertex node in the middle of tet face2
2704 {
2705 unsigned j = 8;
2706
2707 // Need new node?
2709 if (old_node_pt == 0)
2710 {
2711 Node* new_node_pt = 0;
2712 if (face2.is_boundary_face())
2713 {
2715 }
2716 else
2717 {
2719 }
2721 Node_pt.push_back(new_node_pt);
2722 Vector<double> s(3);
2726 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2728 new_node_pt->x(0) = x_tet[0];
2729 new_node_pt->x(1) = x_tet[1];
2730 new_node_pt->x(2) = x_tet[2];
2731 }
2732 // Node already exists
2733 else
2734 {
2736 }
2737 }
2738
2739 // Brick vertex node in centroid of tet. Only built for first element.
2740 // Enumerated "13" in initial sketch.
2741 {
2742 unsigned j = 26;
2743
2744 // Always copied
2746 }
2747
2748 // Internal brick node -- always built
2749 {
2750 unsigned j = 13;
2751
2752 // Always new
2753 {
2755 Node_pt.push_back(new_node_pt);
2756 Vector<double> s(3);
2760 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2762 new_node_pt->x(0) = x_tet[0];
2763 new_node_pt->x(1) = x_tet[1];
2764 new_node_pt->x(2) = x_tet[2];
2765 }
2766 }
2767
2768 // Brick edge node between brick nodes 0 and 2
2769 {
2770 unsigned j = 1;
2771
2772 // Need new node?
2773 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
2775 if (old_node_pt == 0)
2776 {
2777 Node* new_node_pt = 0;
2778 if (edge.is_boundary_edge())
2779 {
2781 }
2782 else
2783 {
2785 }
2787 Node_pt.push_back(new_node_pt);
2788 Vector<double> s(3);
2792 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2794 new_node_pt->x(0) = x_tet[0];
2795 new_node_pt->x(1) = x_tet[1];
2796 new_node_pt->x(2) = x_tet[2];
2797 }
2798 // Node already exists
2799 else
2800 {
2802 }
2803 }
2804
2805 // Brick edge node between brick nodes 0 and 6
2806 {
2807 unsigned j = 3;
2808
2809 // Need new node?
2810 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
2812 if (old_node_pt == 0)
2813 {
2814 Node* new_node_pt = 0;
2815 if (edge.is_boundary_edge())
2816 {
2818 }
2819 else
2820 {
2822 }
2824 Node_pt.push_back(new_node_pt);
2825 Vector<double> s(3);
2829 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2831 new_node_pt->x(0) = x_tet[0];
2832 new_node_pt->x(1) = x_tet[1];
2833 new_node_pt->x(2) = x_tet[2];
2834 }
2835 // Node already exists
2836 else
2837 {
2839 }
2840 }
2841
2842 // Brick edge node between brick nodes 2 and 8
2843 {
2844 unsigned j = 5;
2845
2846 // Need new node?
2847 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
2849 if (old_node_pt == 0)
2850 {
2851 Node* new_node_pt = 0;
2852 if (edge.is_boundary_edge())
2853 {
2855 }
2856 else
2857 {
2859 }
2861 Node_pt.push_back(new_node_pt);
2862 Vector<double> s(3);
2866 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2868 new_node_pt->x(0) = x_tet[0];
2869 new_node_pt->x(1) = x_tet[1];
2870 new_node_pt->x(2) = x_tet[2];
2871 }
2872 // Node already exists
2873 else
2874 {
2876 }
2877 }
2878
2879 // Brick edge node between brick nodes 6 and 8
2880 {
2881 unsigned j = 7;
2882
2883 // Need new node?
2884 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
2886 if (old_node_pt == 0)
2887 {
2888 Node* new_node_pt = 0;
2889 if (edge.is_boundary_edge())
2890 {
2892 }
2893 else
2894 {
2896 }
2898 Node_pt.push_back(new_node_pt);
2899 Vector<double> s(3);
2903 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2905 new_node_pt->x(0) = x_tet[0];
2906 new_node_pt->x(1) = x_tet[1];
2907 new_node_pt->x(2) = x_tet[2];
2908 }
2909 // Node already exists
2910 else
2911 {
2913 }
2914 }
2915
2916 // Brick edge node between brick nodes 18 and 20
2917 {
2918 unsigned j = 19;
2919
2920 // Need new node?
2921 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
2923 if (old_node_pt == 0)
2924 {
2925 Node* new_node_pt = 0;
2926 if (edge.is_boundary_edge())
2927 {
2929 }
2930 else
2931 {
2933 }
2935 Node_pt.push_back(new_node_pt);
2936 Vector<double> s(3);
2940 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2942 new_node_pt->x(0) = x_tet[0];
2943 new_node_pt->x(1) = x_tet[1];
2944 new_node_pt->x(2) = x_tet[2];
2945 }
2946 // Node already exists
2947 else
2948 {
2950 }
2951 }
2952
2953 // Brick edge node between brick nodes 18 and 24
2954 {
2955 unsigned j = 21;
2956
2957 // Need new node?
2958 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
2960 if (old_node_pt == 0)
2961 {
2962 Node* new_node_pt = 0;
2963 if (edge.is_boundary_edge())
2964 {
2966 }
2967 else
2968 {
2970 }
2972 Node_pt.push_back(new_node_pt);
2973 Vector<double> s(3);
2977 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
2979 new_node_pt->x(0) = x_tet[0];
2980 new_node_pt->x(1) = x_tet[1];
2981 new_node_pt->x(2) = x_tet[2];
2982 }
2983 // Node already exists
2984 else
2985 {
2987 }
2988 }
2989
2990 // Brick edge node between brick nodes 20 and 26
2991 {
2992 unsigned j = 23;
2993
2994 // Need new node?
2995 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
2997 if (old_node_pt == 0)
2998 {
2999 Node* new_node_pt = 0;
3000 if (edge.is_boundary_edge())
3001 {
3003 }
3004 else
3005 {
3007 }
3009 Node_pt.push_back(new_node_pt);
3010 Vector<double> s(3);
3014 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3016 new_node_pt->x(0) = x_tet[0];
3017 new_node_pt->x(1) = x_tet[1];
3018 new_node_pt->x(2) = x_tet[2];
3019 }
3020 // Node already exists
3021 else
3022 {
3024 }
3025 }
3026
3027 // Brick edge node between brick nodes 24 and 26
3028 {
3029 unsigned j = 25;
3030
3031 // Need new node?
3032 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
3034 if (old_node_pt == 0)
3035 {
3036 Node* new_node_pt = 0;
3037 if (edge.is_boundary_edge())
3038 {
3040 }
3041 else
3042 {
3044 }
3046 Node_pt.push_back(new_node_pt);
3047 Vector<double> s(3);
3051 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3053 new_node_pt->x(0) = x_tet[0];
3054 new_node_pt->x(1) = x_tet[1];
3055 new_node_pt->x(2) = x_tet[2];
3056 }
3057 // Node already exists
3058 else
3059 {
3061 }
3062 }
3063
3064 // Brick edge node between brick nodes 0 and 18
3065 {
3066 unsigned j = 9;
3067
3068 // Need new node?
3069 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
3071 if (old_node_pt == 0)
3072 {
3073 Node* new_node_pt = 0;
3074 if (edge.is_boundary_edge())
3075 {
3077 }
3078 else
3079 {
3081 }
3083 Node_pt.push_back(new_node_pt);
3084 Vector<double> s(3);
3088 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3090 new_node_pt->x(0) = x_tet[0];
3091 new_node_pt->x(1) = x_tet[1];
3092 new_node_pt->x(2) = x_tet[2];
3093 }
3094 // Node already exists
3095 else
3096 {
3098 }
3099 }
3100
3101 // Brick edge node between brick nodes 2 and 20
3102 {
3103 unsigned j = 11;
3104
3105 // Need new node?
3106 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
3108 if (old_node_pt == 0)
3109 {
3110 Node* new_node_pt = 0;
3111 if (edge.is_boundary_edge())
3112 {
3114 }
3115 else
3116 {
3118 }
3120 Node_pt.push_back(new_node_pt);
3121 Vector<double> s(3);
3125 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3127 new_node_pt->x(0) = x_tet[0];
3128 new_node_pt->x(1) = x_tet[1];
3129 new_node_pt->x(2) = x_tet[2];
3130 }
3131 // Node already exists
3132 else
3133 {
3135 }
3136 }
3137
3138 // Brick edge node between brick nodes 6 and 24
3139 {
3140 unsigned j = 15;
3141
3142 // Need new node?
3143 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
3145 if (old_node_pt == 0)
3146 {
3147 Node* new_node_pt = 0;
3148 if (edge.is_boundary_edge())
3149 {
3151 }
3152 else
3153 {
3155 }
3157 Node_pt.push_back(new_node_pt);
3158 Vector<double> s(3);
3162 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3164 new_node_pt->x(0) = x_tet[0];
3165 new_node_pt->x(1) = x_tet[1];
3166 new_node_pt->x(2) = x_tet[2];
3167 }
3168 // Node already exists
3169 else
3170 {
3172 }
3173 }
3174
3175 // Brick edge node between brick nodes 8 and 26
3176 {
3177 unsigned j = 17;
3178
3179 // Need new node?
3180 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
3182 if (old_node_pt == 0)
3183 {
3184 Node* new_node_pt = 0;
3185 if (edge.is_boundary_edge())
3186 {
3188 }
3189 else
3190 {
3192 }
3194 Node_pt.push_back(new_node_pt);
3195 Vector<double> s(3);
3199 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3201 new_node_pt->x(0) = x_tet[0];
3202 new_node_pt->x(1) = x_tet[1];
3203 new_node_pt->x(2) = x_tet[2];
3204 }
3205 // Node already exists
3206 else
3207 {
3209 }
3210 }
3211
3212 // Mid brick-face node associated with face
3213 // spanned by mid-vertex nodes associated with tet edges 0 and 2
3214 {
3215 unsigned j = 10;
3216
3217 // Need new node?
3221
3223 if (old_node_pt == 0)
3224 {
3225 Node* new_node_pt = 0;
3226 if (face.is_boundary_face())
3227 {
3229 }
3230 else
3231 {
3233 }
3235 Node_pt.push_back(new_node_pt);
3236 Vector<double> s(3);
3240 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3242 new_node_pt->x(0) = x_tet[0];
3243 new_node_pt->x(1) = x_tet[1];
3244 new_node_pt->x(2) = x_tet[2];
3245 }
3246 // Node already exists
3247 else
3248 {
3250 }
3251 }
3252
3253 // Mid brick-face node associated with face
3254 // spanned by mid-vertex nodes associated with tet edges 1 and 2
3255 {
3256 unsigned j = 12;
3257
3258 // Need new node?
3263 if (old_node_pt == 0)
3264 {
3265 Node* new_node_pt = 0;
3266 if (face.is_boundary_face())
3267 {
3269 }
3270 else
3271 {
3273 }
3275 Node_pt.push_back(new_node_pt);
3276 Vector<double> s(3);
3280 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3282 new_node_pt->x(0) = x_tet[0];
3283 new_node_pt->x(1) = x_tet[1];
3284 new_node_pt->x(2) = x_tet[2];
3285 }
3286 // Node already exists
3287 else
3288 {
3290 }
3291 }
3292
3293 // Mid brick-face node associated with face
3294 // spanned by mid-vertex nodes associated with tet edges 0 and 1
3295 {
3296 unsigned j = 4;
3297
3298 // Need new node?
3302
3304 if (old_node_pt == 0)
3305 {
3306 Node* new_node_pt = 0;
3307 if (face.is_boundary_face())
3308 {
3310 }
3311 else
3312 {
3314 }
3316 Node_pt.push_back(new_node_pt);
3317 Vector<double> s(3);
3321 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3323 new_node_pt->x(0) = x_tet[0];
3324 new_node_pt->x(1) = x_tet[1];
3325 new_node_pt->x(2) = x_tet[2];
3326 }
3327 // Node already exists
3328 else
3329 {
3331 }
3332 }
3333
3334 // Top mid brick-face node -- only built by this element
3335 {
3336 unsigned j = 22;
3338 Node_pt.push_back(new_node_pt);
3339 Vector<double> s(3);
3343 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
3346 new_node_pt->x(0) = x_tet[0];
3347 new_node_pt->x(1) = x_tet[1];
3348 new_node_pt->x(2) = x_tet[2];
3349 }
3350
3351 // Right mid brick-face node copy from first element
3352 {
3353 unsigned j = 14;
3354
3355 // Always copied
3357 }
3358
3359 // Back mid brick-face node copy from previous element
3360 {
3361 unsigned j = 16;
3362
3363 // Always copied
3365 }
3366 }
3367
3368 // Fourth brick element is centred at node 2 of tet:
3369 //--------------------------------------------------
3370 {
3371 // Assign coordinates of dummy element
3372 for (unsigned j = 0; j < 8; j++)
3373 {
3377 switch (j)
3378 {
3379 case 0:
3381 nod_pt->set_value(0, s_tet[0]);
3382 nod_pt->set_value(1, s_tet[1]);
3383 nod_pt->set_value(2, s_tet[2]);
3385 nod_pt->x(0) = x_tet[0];
3386 nod_pt->x(1) = x_tet[1];
3387 nod_pt->x(2) = x_tet[2];
3388 break;
3389 case 1:
3391 nod_pt->set_value(0, s_tet[0]);
3392 nod_pt->set_value(1, s_tet[1]);
3393 nod_pt->set_value(2, s_tet[2]);
3395 nod_pt->x(0) = x_tet[0];
3396 nod_pt->x(1) = x_tet[1];
3397 nod_pt->x(2) = x_tet[2];
3398 break;
3399 case 2:
3401 nod_pt->set_value(0, s_tet[0]);
3402 nod_pt->set_value(1, s_tet[1]);
3403 nod_pt->set_value(2, s_tet[2]);
3405 nod_pt->x(0) = x_tet[0];
3406 nod_pt->x(1) = x_tet[1];
3407 nod_pt->x(2) = x_tet[2];
3408 break;
3409 case 3:
3410 // label 11 in initial sketch: Mid face node on face
3411 // spanned by tet nodes 0,1,2
3412 s_tet[0] = 1.0 / 3.0;
3413 s_tet[1] = 1.0 / 3.0;
3414 s_tet[2] = 1.0 / 3.0;
3415 nod_pt->set_value(0, s_tet[0]);
3416 nod_pt->set_value(1, s_tet[1]);
3417 nod_pt->set_value(2, s_tet[2]);
3419 nod_pt->x(0) = x_tet[0];
3420 nod_pt->x(1) = x_tet[1];
3421 nod_pt->x(2) = x_tet[2];
3422 break;
3423 case 4:
3425 nod_pt->set_value(0, s_tet[0]);
3426 nod_pt->set_value(1, s_tet[1]);
3427 nod_pt->set_value(2, s_tet[2]);
3429 nod_pt->x(0) = x_tet[0];
3430 nod_pt->x(1) = x_tet[1];
3431 nod_pt->x(2) = x_tet[2];
3432 break;
3433 case 5:
3434 // label 10 in initial sketch: Mid face node on face
3435 // spanned by tet nodes 0,2,3
3436 s_tet[0] = 0.0;
3437 s_tet[1] = 1.0 / 3.0;
3438 s_tet[2] = 1.0 / 3.0;
3439 nod_pt->set_value(0, s_tet[0]);
3440 nod_pt->set_value(1, s_tet[1]);
3441 nod_pt->set_value(2, s_tet[2]);
3443 nod_pt->x(0) = x_tet[0];
3444 nod_pt->x(1) = x_tet[1];
3445 nod_pt->x(2) = x_tet[2];
3446 break;
3447 case 6:
3448 // label 12 in initial sketch: Mid face node on face
3449 // spanned by tet nodes 0,2,3
3450 s_tet[0] = 1.0 / 3.0;
3451 s_tet[1] = 0.0;
3452 s_tet[2] = 1.0 / 3.0;
3453 nod_pt->set_value(0, s_tet[0]);
3454 nod_pt->set_value(1, s_tet[1]);
3455 nod_pt->set_value(2, s_tet[2]);
3457 nod_pt->x(0) = x_tet[0];
3458 nod_pt->x(1) = x_tet[1];
3459 nod_pt->x(2) = x_tet[2];
3460 break;
3461 case 7:
3462 // label 14 in initial sketch: Centroid
3463 s_tet[0] = 0.25;
3464 s_tet[1] = 0.25;
3465 s_tet[2] = 0.25;
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]);
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 }
3475 }
3476
3477 // Create actual third brick element
3478 FiniteElement* el_pt = new ELEMENT;
3480 Element_pt.push_back(el_pt);
3481
3482 TFace face0(
3484
3485 TFace face1(
3487
3488 TFace face2(
3490
3491 // Tet vertex nodes along edges emanating from node 0 in brick
3493 tet_edge_node[0].resize(2);
3494 tet_edge_node[0][0] = 7;
3495 tet_edge_node[0][1] = 1;
3496 tet_edge_node[1].resize(2);
3497 tet_edge_node[1][0] = 5;
3498 tet_edge_node[1][1] = 0;
3499 tet_edge_node[2].resize(2);
3500 tet_edge_node[2][0] = 8;
3501 tet_edge_node[2][1] = 3;
3502
3503 // Node number of tet vertex that node 0 in brick is centred on
3504 unsigned central_tet_vertex = 2;
3505
3506 Node* tet_node_pt = 0;
3507 Node* old_node_pt = 0;
3508
3509 // Corner node
3510 {
3511 unsigned j = 0;
3512
3513 // Need new node?
3516 if (old_node_pt == 0)
3517 {
3518 Node* new_node_pt = 0;
3519 if (tet_node_pt->is_on_boundary())
3520 {
3522 }
3523 else
3524 {
3526 }
3528 Node_pt.push_back(new_node_pt);
3529 Vector<double> s(3);
3533 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3535 new_node_pt->x(0) = x_tet[0];
3536 new_node_pt->x(1) = x_tet[1];
3537 new_node_pt->x(2) = x_tet[2];
3538 }
3539 // Node already exists
3540 else
3541 {
3543 }
3544 }
3545
3546 // Brick vertex node coindides with mid-edge node on tet edge 0
3547 {
3548 unsigned j = 2;
3549
3550 // Need new node?
3553 if (old_node_pt == 0)
3554 {
3555 Node* new_node_pt = 0;
3556 if (tet_node_pt->is_on_boundary())
3557 {
3559 }
3560 else
3561 {
3563 }
3565 Node_pt.push_back(new_node_pt);
3566 Vector<double> s(3);
3570 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3572 new_node_pt->x(0) = x_tet[0];
3573 new_node_pt->x(1) = x_tet[1];
3574 new_node_pt->x(2) = x_tet[2];
3575 }
3576 // Node already exists
3577 else
3578 {
3580 }
3581 }
3582
3583 // Brick vertex node coindides with mid vertex node of tet edge 1
3584 {
3585 unsigned j = 6;
3586
3587 // Need new node?
3590 if (old_node_pt == 0)
3591 {
3592 Node* new_node_pt = 0;
3593 if (tet_node_pt->is_on_boundary())
3594 {
3596 }
3597 else
3598 {
3600 }
3602 Node_pt.push_back(new_node_pt);
3603 Vector<double> s(3);
3607 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3609 new_node_pt->x(0) = x_tet[0];
3610 new_node_pt->x(1) = x_tet[1];
3611 new_node_pt->x(2) = x_tet[2];
3612 }
3613 // Node already exists
3614 else
3615 {
3617 }
3618 }
3619
3620 // Brick vertex node coindides with mid-vertex node of tet edge 2
3621 {
3622 unsigned j = 18;
3623
3624 // Need new node?
3627 if (old_node_pt == 0)
3628 {
3629 Node* new_node_pt = 0;
3630 if (tet_node_pt->is_on_boundary())
3631 {
3633 }
3634 else
3635 {
3637 }
3639 Node_pt.push_back(new_node_pt);
3640 Vector<double> s(3);
3644 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3646 new_node_pt->x(0) = x_tet[0];
3647 new_node_pt->x(1) = x_tet[1];
3648 new_node_pt->x(2) = x_tet[2];
3649 }
3650 // Node already exists
3651 else
3652 {
3654 }
3655 }
3656
3657 // Brick vertex node in the middle of tet face0
3658 {
3659 unsigned j = 20;
3660
3661 // Need new node?
3663 if (old_node_pt == 0)
3664 {
3665 Node* new_node_pt = 0;
3666 if (face0.is_boundary_face())
3667 {
3669 }
3670 else
3671 {
3673 }
3675 Node_pt.push_back(new_node_pt);
3676 Vector<double> s(3);
3680 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3682 new_node_pt->x(0) = x_tet[0];
3683 new_node_pt->x(1) = x_tet[1];
3684 new_node_pt->x(2) = x_tet[2];
3685 }
3686 // Node already exists
3687 else
3688 {
3690 }
3691 }
3692
3693 // Brick vertex node in the middle of tet face1
3694 {
3695 unsigned j = 24;
3696
3697 // Need new node?
3699 if (old_node_pt == 0)
3700 {
3701 Node* new_node_pt = 0;
3702 if (face1.is_boundary_face())
3703 {
3705 }
3706 else
3707 {
3709 }
3711 Node_pt.push_back(new_node_pt);
3712 Vector<double> s(3);
3716 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3718 new_node_pt->x(0) = x_tet[0];
3719 new_node_pt->x(1) = x_tet[1];
3720 new_node_pt->x(2) = x_tet[2];
3721 }
3722 // Node already exists
3723 else
3724 {
3726 }
3727 }
3728
3729 // Brick vertex node in the middle of tet face2
3730 {
3731 unsigned j = 8;
3732
3733 // Need new node?
3735 if (old_node_pt == 0)
3736 {
3737 Node* new_node_pt = 0;
3738 if (face2.is_boundary_face())
3739 {
3741 }
3742 else
3743 {
3745 }
3747 Node_pt.push_back(new_node_pt);
3748 Vector<double> s(3);
3752 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3754 new_node_pt->x(0) = x_tet[0];
3755 new_node_pt->x(1) = x_tet[1];
3756 new_node_pt->x(2) = x_tet[2];
3757 }
3758 // Node already exists
3759 else
3760 {
3762 }
3763 }
3764
3765 // Brick vertex node in centroid of tet. Only built for first element.
3766 // Enumerated "13" in initial sketch.
3767 {
3768 unsigned j = 26;
3769
3770 // Always copied
3772 }
3773
3774 // Internal brick node -- always built
3775 {
3776 unsigned j = 13;
3777
3778 // Always new
3779 {
3781 Node_pt.push_back(new_node_pt);
3782 Vector<double> s(3);
3786 dummy_q_el_pt[3]->interpolated_s_tet(s, s_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 }
3793
3794 // Brick edge node between brick nodes 0 and 2
3795 {
3796 unsigned j = 1;
3797
3798 // Need new node?
3799 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
3801 if (old_node_pt == 0)
3802 {
3803 Node* new_node_pt = 0;
3804 if (edge.is_boundary_edge())
3805 {
3807 }
3808 else
3809 {
3811 }
3813 Node_pt.push_back(new_node_pt);
3814 Vector<double> s(3);
3818 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3820 new_node_pt->x(0) = x_tet[0];
3821 new_node_pt->x(1) = x_tet[1];
3822 new_node_pt->x(2) = x_tet[2];
3823 }
3824 // Node already exists
3825 else
3826 {
3828 }
3829 }
3830
3831 // Brick edge node between brick nodes 0 and 6
3832 {
3833 unsigned j = 3;
3834
3835 // Need new node?
3836 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
3838 if (old_node_pt == 0)
3839 {
3840 Node* new_node_pt = 0;
3841 if (edge.is_boundary_edge())
3842 {
3844 }
3845 else
3846 {
3848 }
3850 Node_pt.push_back(new_node_pt);
3851 Vector<double> s(3);
3855 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3857 new_node_pt->x(0) = x_tet[0];
3858 new_node_pt->x(1) = x_tet[1];
3859 new_node_pt->x(2) = x_tet[2];
3860 }
3861 // Node already exists
3862 else
3863 {
3865 }
3866 }
3867
3868 // Brick edge node between brick nodes 2 and 8
3869 {
3870 unsigned j = 5;
3871
3872 // Need new node?
3873 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
3875 if (old_node_pt == 0)
3876 {
3877 Node* new_node_pt = 0;
3878 if (edge.is_boundary_edge())
3879 {
3881 }
3882 else
3883 {
3885 }
3887 Node_pt.push_back(new_node_pt);
3888 Vector<double> s(3);
3892 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3894 new_node_pt->x(0) = x_tet[0];
3895 new_node_pt->x(1) = x_tet[1];
3896 new_node_pt->x(2) = x_tet[2];
3897 }
3898 // Node already exists
3899 else
3900 {
3902 }
3903 }
3904
3905 // Brick edge node between brick nodes 6 and 8
3906 {
3907 unsigned j = 7;
3908
3909 // Need new node?
3910 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
3912 if (old_node_pt == 0)
3913 {
3914 Node* new_node_pt = 0;
3915 if (edge.is_boundary_edge())
3916 {
3918 }
3919 else
3920 {
3922 }
3924 Node_pt.push_back(new_node_pt);
3925 Vector<double> s(3);
3929 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3931 new_node_pt->x(0) = x_tet[0];
3932 new_node_pt->x(1) = x_tet[1];
3933 new_node_pt->x(2) = x_tet[2];
3934 }
3935 // Node already exists
3936 else
3937 {
3939 }
3940 }
3941
3942 // Brick edge node between brick nodes 18 and 20
3943 {
3944 unsigned j = 19;
3945
3946 // Need new node?
3947 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
3949 if (old_node_pt == 0)
3950 {
3951 Node* new_node_pt = 0;
3952 if (edge.is_boundary_edge())
3953 {
3955 }
3956 else
3957 {
3959 }
3961 Node_pt.push_back(new_node_pt);
3962 Vector<double> s(3);
3966 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
3968 new_node_pt->x(0) = x_tet[0];
3969 new_node_pt->x(1) = x_tet[1];
3970 new_node_pt->x(2) = x_tet[2];
3971 }
3972 // Node already exists
3973 else
3974 {
3976 }
3977 }
3978
3979 // Brick edge node between brick nodes 18 and 24
3980 {
3981 unsigned j = 21;
3982
3983 // Need new node?
3984 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
3986 if (old_node_pt == 0)
3987 {
3988 Node* new_node_pt = 0;
3989 if (edge.is_boundary_edge())
3990 {
3992 }
3993 else
3994 {
3996 }
3998 Node_pt.push_back(new_node_pt);
3999 Vector<double> s(3);
4003 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
4005 new_node_pt->x(0) = x_tet[0];
4006 new_node_pt->x(1) = x_tet[1];
4007 new_node_pt->x(2) = x_tet[2];
4008 }
4009 // Node already exists
4010 else
4011 {
4013 }
4014 }
4015
4016 // Brick edge node between brick nodes 20 and 26
4017 {
4018 unsigned j = 23;
4019
4020 // Need new node?
4021 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
4023 if (old_node_pt == 0)
4024 {
4025 Node* new_node_pt = 0;
4026 if (edge.is_boundary_edge())
4027 {
4029 }
4030 else
4031 {
4033 }
4035 Node_pt.push_back(new_node_pt);
4036 Vector<double> s(3);
4040 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
4042 new_node_pt->x(0) = x_tet[0];
4043 new_node_pt->x(1) = x_tet[1];
4044 new_node_pt->x(2) = x_tet[2];
4045 }
4046 // Node already exists
4047 else
4048 {
4050 }
4051 }
4052
4053 // Brick edge node between brick nodes 24 and 26
4054 {
4055 unsigned j = 25;
4056
4057 // Need new node?
4058 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
4060 if (old_node_pt == 0)
4061 {
4062 Node* new_node_pt = 0;
4063 if (edge.is_boundary_edge())
4064 {
4066 }
4067 else
4068 {
4070 }
4072 Node_pt.push_back(new_node_pt);
4073 Vector<double> s(3);
4077 dummy_q_el_pt[3]->interpolated_s_tet(s, s_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 {
4087 }
4088 }
4089
4090 // Brick edge node between brick nodes 0 and 18
4091 {
4092 unsigned j = 9;
4093
4094 // Need new node?
4095 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
4097 if (old_node_pt == 0)
4098 {
4099 Node* new_node_pt = 0;
4100 if (edge.is_boundary_edge())
4101 {
4103 }
4104 else
4105 {
4107 }
4109 Node_pt.push_back(new_node_pt);
4110 Vector<double> s(3);
4114 dummy_q_el_pt[3]->interpolated_s_tet(s, s_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 {
4124 }
4125 }
4126
4127 // Brick edge node between brick nodes 2 and 20
4128 {
4129 unsigned j = 11;
4130
4131 // Need new node?
4132 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
4134 if (old_node_pt == 0)
4135 {
4136 Node* new_node_pt = 0;
4137 if (edge.is_boundary_edge())
4138 {
4140 }
4141 else
4142 {
4144 }
4146 Node_pt.push_back(new_node_pt);
4147 Vector<double> s(3);
4151 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
4153 new_node_pt->x(0) = x_tet[0];
4154 new_node_pt->x(1) = x_tet[1];
4155 new_node_pt->x(2) = x_tet[2];
4156 }
4157 // Node already exists
4158 else
4159 {
4161 }
4162 }
4163
4164 // Brick edge node between brick nodes 6 and 24
4165 {
4166 unsigned j = 15;
4167
4168 // Need new node?
4169 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
4171 if (old_node_pt == 0)
4172 {
4173 Node* new_node_pt = 0;
4174 if (edge.is_boundary_edge())
4175 {
4177 }
4178 else
4179 {
4181 }
4183 Node_pt.push_back(new_node_pt);
4184 Vector<double> s(3);
4188 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
4190 new_node_pt->x(0) = x_tet[0];
4191 new_node_pt->x(1) = x_tet[1];
4192 new_node_pt->x(2) = x_tet[2];
4193 }
4194 // Node already exists
4195 else
4196 {
4198 }
4199 }
4200
4201 // Brick edge node between brick nodes 8 and 26
4202 {
4203 unsigned j = 17;
4204
4205 // Need new node?
4206 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
4208 if (old_node_pt == 0)
4209 {
4210 Node* new_node_pt = 0;
4211 if (edge.is_boundary_edge())
4212 {
4214 }
4215 else
4216 {
4218 }
4220 Node_pt.push_back(new_node_pt);
4221 Vector<double> s(3);
4225 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
4227 new_node_pt->x(0) = x_tet[0];
4228 new_node_pt->x(1) = x_tet[1];
4229 new_node_pt->x(2) = x_tet[2];
4230 }
4231 // Node already exists
4232 else
4233 {
4235 }
4236 }
4237
4238 // Mid brick-face node associated with face
4239 // spanned by mid-vertex nodes associated with tet edges 0 and 2
4240 {
4241 unsigned j = 10;
4242
4243 // Need new node?
4247
4249 if (old_node_pt == 0)
4250 {
4251 Node* new_node_pt = 0;
4252 if (face.is_boundary_face())
4253 {
4255 }
4256 else
4257 {
4259 }
4261 Node_pt.push_back(new_node_pt);
4262 Vector<double> s(3);
4266 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
4268 new_node_pt->x(0) = x_tet[0];
4269 new_node_pt->x(1) = x_tet[1];
4270 new_node_pt->x(2) = x_tet[2];
4271 }
4272 // Node already exists
4273 else
4274 {
4276 }
4277 }
4278
4279 // Mid brick-face node associated with face
4280 // spanned by mid-vertex nodes associated with tet edges 1 and 2
4281 {
4282 unsigned j = 12;
4283
4284 // Need new node?
4288
4290 if (old_node_pt == 0)
4291 {
4292 Node* new_node_pt = 0;
4293 if (face.is_boundary_face())
4294 {
4296 }
4297 else
4298 {
4300 }
4302 Node_pt.push_back(new_node_pt);
4303 Vector<double> s(3);
4307 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
4309 new_node_pt->x(0) = x_tet[0];
4310 new_node_pt->x(1) = x_tet[1];
4311 new_node_pt->x(2) = x_tet[2];
4312 }
4313 // Node already exists
4314 else
4315 {
4317 }
4318 }
4319
4320 // Mid brick-face node associated with face
4321 // spanned by mid-vertex nodes associated with tet edges 0 and 1
4322 {
4323 unsigned j = 4;
4324
4325 // Need new node?
4329
4331 if (old_node_pt == 0)
4332 {
4333 Node* new_node_pt = 0;
4334 if (face.is_boundary_face())
4335 {
4337 }
4338 else
4339 {
4341 }
4343 Node_pt.push_back(new_node_pt);
4344 Vector<double> s(3);
4348 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
4350 new_node_pt->x(0) = x_tet[0];
4351 new_node_pt->x(1) = x_tet[1];
4352 new_node_pt->x(2) = x_tet[2];
4353 }
4354 // Node already exists
4355 else
4356 {
4358 }
4359 }
4360
4361 // Top mid brick-face node copy from top of second element
4362 {
4363 unsigned j = 22;
4364
4365 // Always copied
4367 }
4368
4369 // Right mid brick-face node copy from top of first element
4370 {
4371 unsigned j = 14;
4372
4373 // Always copied
4375 }
4376
4377 // Back mid brick-face node copy from top of zeroth element
4378 {
4379 unsigned j = 16;
4380
4381 // Always copied
4383 }
4384 }
4385
4386 // Check if the four faces of the tet are on a boundary.
4387 // If they are, add the nodes in the brick mesh to those
4388 // boundaries.
4389
4390 // Loop over faces of tet (in its own face index enumeration)
4391 for (int face_index = 0; face_index < 4; face_index++)
4392 {
4393 // Identify face and coordinates in it
4394 TFace* face_pt = 0;
4395 switch (face_index)
4396 {
4397 case 0:
4398 // Face 0: s[0]=0; opposite node 0
4399 face_pt = new TFace(tet_el_pt->node_pt(1),
4400 tet_el_pt->node_pt(2),
4401 tet_el_pt->node_pt(3));
4402 break;
4403
4404 case 1:
4405 // Face 1: s[1]=0; opposite node 1
4406 face_pt = new TFace(tet_el_pt->node_pt(0),
4407 tet_el_pt->node_pt(2),
4408 tet_el_pt->node_pt(3));
4409
4410 break;
4411
4412 case 2:
4413 // Face 2: s[2]=0; opposite node 2
4414 face_pt = new TFace(tet_el_pt->node_pt(0),
4415 tet_el_pt->node_pt(1),
4416 tet_el_pt->node_pt(3));
4417
4418 break;
4419
4420 case 3:
4421 // Face 3: s[0]+s[1]+s[2]=1; opposite node 3
4422 face_pt = new TFace(tet_el_pt->node_pt(0),
4423 tet_el_pt->node_pt(1),
4424 tet_el_pt->node_pt(2));
4425 break;
4426 }
4427
4428 if (face_pt->is_boundary_face())
4429 {
4430 std::set<unsigned> bnd;
4431 std::set<unsigned>* bnd_pt = &bnd;
4432 face_pt->get_boundaries_pt(bnd_pt);
4433
4434#ifdef PARANOID
4435 if ((*bnd_pt).size() > 1)
4436 {
4437 std::ostringstream error_stream;
4438 error_stream << "TFace should only be on one boundary.\n";
4439 throw OomphLibError(error_stream.str(),
4442 }
4443#endif
4444
4445 if ((*bnd_pt).size() == 1)
4446 {
4447 // Create new face element
4450 {
4451#ifdef PARANOID
4452 if (dynamic_cast<SolidTElement<3, 3>*>(tet_el_pt) == 0)
4453 {
4454 std::ostringstream error_stream;
4456 << "Tet-element cannot be cast to SolidTElement<3,3>.\n"
4457 << "BrickFromTetMesh can only be built from\n"
4458 << "mesh containing quadratic tets.\n"
4459 << std::endl;
4460 throw OomphLibError(error_stream.str(),
4463 }
4464#endif
4465 // Build the face element
4467 tet_el_pt, face_index);
4468 }
4469 else
4470 {
4471#ifdef PARANOID
4472 if (dynamic_cast<TElement<3, 3>*>(tet_el_pt) == 0)
4473 {
4474 std::ostringstream error_stream;
4475 error_stream << "Tet-element cannot be cast to TElement<3,3>.\n"
4476 << "BrickFromTetMesh can only be built from\n"
4477 << "mesh containing quadratic tets.\n"
4478 << std::endl;
4479 throw OomphLibError(error_stream.str(),
4482 }
4483#endif
4484
4485 // Build the face element
4486 face_el_pt =
4488 }
4489
4490 // Specify boundary id in bulk mesh (needed to extract
4491 // boundary coordinate)
4492 unsigned b = (*(*bnd_pt).begin());
4493 set_boundary_coordinate_exists(b);
4494 face_el_pt->set_boundary_number_in_bulk_mesh(b);
4495
4496 // Now set up the brick nodes on this face, enumerated as
4497 // in s_face
4499
4500 switch (face_index)
4501 {
4502 case 0:
4506
4510
4514
4518
4520
4523
4526
4529 break;
4530
4531 case 1:
4535
4539
4543
4547
4549
4552
4555
4558 break;
4559
4560 case 2:
4564
4568
4572
4576
4578
4581
4584
4587 break;
4588
4589 case 3:
4593
4597
4601
4605
4607
4610
4613
4616 break;
4617 }
4618
4619 // Provide possibility for translation -- may need to add
4620 // face index to this; see ThinLayerBrickOnTetMesh.
4622
4623 // Initialise with identity mapping
4624 for (unsigned i = 0; i < 19; i++)
4625 {
4626 translate[i] = i;
4627 }
4628
4629 // Visit all the nodes on that face
4630 for (unsigned j = 0; j < 19; j++)
4631 {
4632 // Which node is it?
4634
4635 // Get coordinates etc of point from face
4638 Vector<double> x(3);
4641
4642#ifdef PARANOID
4643 // Check that the coordinates match (within tolerance)
4644 double dist = sqrt(pow(brick_node_pt->x(0) - x[0], 2) +
4645 pow(brick_node_pt->x(1) - x[1], 2) +
4646 pow(brick_node_pt->x(2) - x[2], 2));
4648 {
4649 std::ofstream brick0;
4650 std::ofstream brick1;
4651 std::ofstream brick2;
4652 std::ofstream brick3;
4653 brick0.open("full_brick0.dat");
4654 brick1.open("full_brick1.dat");
4655 brick2.open("full_brick2.dat");
4656 brick3.open("full_brick3.dat");
4657 for (unsigned j = 0; j < 27; j++)
4658 {
4659 brick0 << brick_el0_pt->node_pt(j)->x(0) << " "
4660 << brick_el0_pt->node_pt(j)->x(1) << " "
4661 << brick_el0_pt->node_pt(j)->x(2) << "\n";
4662
4663 brick1 << brick_el1_pt->node_pt(j)->x(0) << " "
4664 << brick_el1_pt->node_pt(j)->x(1) << " "
4665 << brick_el1_pt->node_pt(j)->x(2) << "\n";
4666
4667 brick2 << brick_el2_pt->node_pt(j)->x(0) << " "
4668 << brick_el2_pt->node_pt(j)->x(1) << " "
4669 << brick_el2_pt->node_pt(j)->x(2) << "\n";
4670
4671 brick3 << brick_el3_pt->node_pt(j)->x(0) << " "
4672 << brick_el3_pt->node_pt(j)->x(1) << " "
4673 << brick_el3_pt->node_pt(j)->x(2) << "\n";
4674 }
4675 brick0.close();
4676 brick1.close();
4677 brick2.close();
4678 brick3.close();
4679
4680 std::ofstream full_face;
4681 full_face.open("full_face.dat");
4682 for (unsigned j = 0; j < 6; j++)
4683 {
4684 full_face << face_el_pt->node_pt(j)->x(0) << " "
4685 << face_el_pt->node_pt(j)->x(1) << " "
4686 << face_el_pt->node_pt(j)->x(2) << "\n";
4687 }
4688 full_face.close();
4689
4690 // Get normal sign
4691 int normal_sign = face_el_pt->normal_sign();
4692
4693 std::ostringstream error_stream;
4695 << "During assignment of boundary cordinates, the distance\n"
4696 << "between brick node and reference point in \n"
4697 << "triangular FaceElement is " << dist << " which \n"
4698 << "is bigger than the tolerance defined in \n"
4699 << "BrickFromTetMeshHelper::Face_position_tolerance="
4701 << "If this is tolerable, increase the tolerance \n"
4702 << "(it's defined in a namespace and therefore publically\n"
4703 << "accessible). If not, the Face may be inverted in which \n"
4704 << "case you should re-implement the translation scheme,\n"
4705 << "following the pattern used in the "
4706 "ThinLayerBrickOnTetMesh."
4707 << "\nThe required code fragements are already here but \n"
4708 << "the translation is the unit map.\n"
4709 << "To aid the diagnostics, the files full_brick[0-3].dat\n"
4710 << "contain the coordinates of the 27 nodes in the four\n"
4711 << "bricks associated with the current tet and "
4712 "full_face.dat\n"
4713 << "contains the coordinates of the 6 nodes in the "
4714 "FaceElement"
4715 << "\nfrom which the boundary coordinates are extracted.\n"
4716 << "FYI: The normal_sign of the face is: " << normal_sign
4717 << std::endl;
4718 throw OomphLibError(error_stream.str(),
4721 }
4722#endif
4723
4724 // Set boundary stuff
4725 add_boundary_node(b, brick_node_pt);
4726 brick_node_pt->set_coordinates_on_boundary(b, zeta);
4727 }
4728
4729 // Add appropriate brick elements to boundary lookup scheme
4730 switch (face_index)
4731 {
4732 case 0:
4733 Boundary_element_pt[b].push_back(brick_el1_pt);
4734 Face_index_at_boundary[b].push_back(-2);
4735 Boundary_element_pt[b].push_back(brick_el2_pt);
4736 Face_index_at_boundary[b].push_back(-1);
4737 Boundary_element_pt[b].push_back(brick_el3_pt);
4738 Face_index_at_boundary[b].push_back(-2);
4739 break;
4740
4741 case 1:
4742 Boundary_element_pt[b].push_back(brick_el0_pt);
4743 Face_index_at_boundary[b].push_back(-1);
4744 Boundary_element_pt[b].push_back(brick_el2_pt);
4745 Face_index_at_boundary[b].push_back(-2);
4746 Boundary_element_pt[b].push_back(brick_el3_pt);
4747 Face_index_at_boundary[b].push_back(-1);
4748 break;
4749
4750 case 2:
4751 Boundary_element_pt[b].push_back(brick_el0_pt);
4752 Face_index_at_boundary[b].push_back(-3);
4753 Boundary_element_pt[b].push_back(brick_el1_pt);
4754 Face_index_at_boundary[b].push_back(-3);
4755 Boundary_element_pt[b].push_back(brick_el2_pt);
4756 Face_index_at_boundary[b].push_back(-3);
4757 break;
4758
4759 case 3:
4760 Boundary_element_pt[b].push_back(brick_el0_pt);
4761 Face_index_at_boundary[b].push_back(-2);
4762 Boundary_element_pt[b].push_back(brick_el1_pt);
4763 Face_index_at_boundary[b].push_back(-1);
4764 Boundary_element_pt[b].push_back(brick_el3_pt);
4765 Face_index_at_boundary[b].push_back(-3);
4766 break;
4767 }
4768 // Cleanup
4769 delete face_el_pt;
4770 }
4771 }
4772 // Cleanup
4773 delete face_pt;
4774 }
4775 }
4776
4777 // Lookup scheme has now been setup
4778 Lookup_for_elements_next_boundary_is_setup = true;
4779
4780 // Get number of distinct boundaries specified
4781 // in the original xda enumeration.
4782 unsigned n_xda_boundaries = tet_mesh_pt->nxda_boundary();
4783
4784 // Copy collective IDs across
4785 Boundary_id.resize(n_xda_boundaries);
4786 for (unsigned xda_b = 0; xda_b < n_xda_boundaries; xda_b++)
4787 {
4788 Boundary_id[xda_b] = tet_mesh_pt->oomph_lib_boundary_ids(xda_b);
4789 }
4790
4791 // Cleanup
4792 for (unsigned e = 0; e < 4; e++)
4793 {
4794 for (unsigned j = 0; j < 8; j++)
4795 {
4796 delete dummy_q_el_pt[e]->node_pt(j);
4797 }
4798 delete dummy_q_el_pt[e];
4799 }
4800 }
4801
4802 //=======================================================================
4803 /// Build fct: Pass pointer to existing tet mesh and timestepper
4804 /// Specialisation for TetgenMesh<TElement<3,3> >
4805 //=======================================================================
4806 template<class ELEMENT>
4809 {
4810 // Mesh can only be built with 3D Qelements.
4811 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3, 3);
4812
4813 // Figure out if the tet mesh is a solid mesh
4814 bool tet_mesh_is_solid_mesh = false;
4815 if (dynamic_cast<SolidFiniteElement*>(tet_mesh_pt->element_pt(0)) != 0)
4816 {
4818 }
4819
4820 // Setup lookup scheme for local coordinates on triangular faces.
4821 // The local coordinates identify the points on the triangular
4822 // FaceElements on which we place the bottom layer of the
4823 // brick nodes.
4825 for (unsigned i = 0; i < 19; i++)
4826 {
4827 s_face[i].resize(2);
4828
4829 switch (i)
4830 {
4831 // Vertex nodes
4832
4833 case 0:
4834 s_face[i][0] = 1.0;
4835 s_face[i][1] = 0.0;
4836 break;
4837
4838 case 1:
4839 s_face[i][0] = 0.0;
4840 s_face[i][1] = 1.0;
4841 break;
4842
4843 case 2:
4844 s_face[i][0] = 0.0;
4845 s_face[i][1] = 0.0;
4846 break;
4847
4848 // Midside nodes
4849
4850 case 3:
4851 s_face[i][0] = 0.5;
4852 s_face[i][1] = 0.5;
4853 break;
4854
4855 case 4:
4856 s_face[i][0] = 0.0;
4857 s_face[i][1] = 0.5;
4858 break;
4859
4860 case 5:
4861 s_face[i][0] = 0.5;
4862 s_face[i][1] = 0.0;
4863 break;
4864
4865 // Quarter side nodes
4866
4867 case 6:
4868 s_face[i][0] = 0.75;
4869 s_face[i][1] = 0.25;
4870 break;
4871
4872 case 7:
4873 s_face[i][0] = 0.25;
4874 s_face[i][1] = 0.75;
4875 break;
4876
4877 case 8:
4878 s_face[i][0] = 0.0;
4879 s_face[i][1] = 0.75;
4880 break;
4881
4882 case 9:
4883 s_face[i][0] = 0.0;
4884 s_face[i][1] = 0.25;
4885 break;
4886
4887 case 10:
4888 s_face[i][0] = 0.25;
4889 s_face[i][1] = 0.0;
4890 break;
4891
4892 case 11:
4893 s_face[i][0] = 0.75;
4894 s_face[i][1] = 0.0;
4895 break;
4896
4897 // Central node
4898
4899 case 12:
4900 s_face[i][0] = 1.0 / 3.0;
4901 s_face[i][1] = 1.0 / 3.0;
4902 break;
4903
4904 // Vertical internal midside nodes connecting 2 and 3
4905
4906 case 13:
4907 s_face[i][0] = 5.0 / 24.0;
4908 s_face[i][1] = 5.0 / 24.0;
4909 break;
4910
4911 case 14:
4912 s_face[i][0] = 5.0 / 12.0;
4913 s_face[i][1] = 5.0 / 12.0;
4914 break;
4915
4916 // Internal midside nodes connecting nodes 0 and 4
4917
4918 case 15:
4919 s_face[i][1] = 5.0 / 24.0;
4920 s_face[i][0] = 7.0 / 12.0; // 1.0-2.0*5.0/24.0;
4921 break;
4922
4923 case 16:
4924 s_face[i][1] = 5.0 / 12.0;
4925 s_face[i][0] = 1.0 / 6.0; // 1.0-2.0*5.0/12.0;
4926 break;
4927
4928 // Internal midside nodes connecting nodes 1 and 5
4929
4930 case 17:
4931 s_face[i][0] = 5.0 / 24.0;
4932 s_face[i][1] = 7.0 / 12.0; // 1.0-2.0*5.0/24.0;
4933 break;
4934
4935 case 18:
4936 s_face[i][0] = 5.0 / 12.0;
4937 s_face[i][1] = 1.0 / 6.0; // 1.0-2.0*5.0/12.0;
4938 break;
4939 }
4940 }
4941
4942 // Set number of boundaries
4943 unsigned nb = tet_mesh_pt->nboundary();
4944 set_nboundary(nb);
4945
4946 // Get ready for boundary lookup scheme
4947 Boundary_element_pt.resize(nb);
4948 Face_index_at_boundary.resize(nb);
4949
4950 // Maps to check which nodes have already been done
4951
4952 // Map that stores the new brick node corresponding to an existing tet node
4953 std::map<Node*, Node*> tet_node_node_pt;
4954
4955 // Map that stores node on an edge between two brick nodes
4956 std::map<Edge, Node*> brick_edge_node_pt;
4957
4958 // Map that stores node on face spanned by three tet nodes
4959 std::map<TFace, Node*> tet_face_node_pt;
4960
4961 // Create the four Dummy bricks:
4962 //------------------------------
4964 for (unsigned e = 0; e < 4; e++)
4965 {
4967 for (unsigned j = 0; j < 8; j++)
4968 {
4970 }
4971 }
4972
4973 // Loop over the elements in the tet mesh
4974 unsigned n_el_tet = tet_mesh_pt->nelement();
4975 for (unsigned e_tet = 0; e_tet < n_el_tet; e_tet++)
4976 {
4977 // Cast to ten-noded tet
4979 dynamic_cast<TElement<3, 3>*>(tet_mesh_pt->element_pt(e_tet));
4980
4981#ifdef PARANOID
4982 if (tet_el_pt == 0)
4983 {
4984 std::ostringstream error_stream;
4986 << "BrickFromTetMesh can only built from tet mesh containing\n"
4987 << "ten-noded tets.\n";
4988 throw OomphLibError(
4990 }
4991#endif
4992
4993 // Storage for the centroid node for this tet
4995
4996 // Internal mid brick-face nodes
5000
5003
5005
5006 // Newly created brick elements
5011
5012 // First brick element is centred at node 0 of tet:
5013 //-------------------------------------------------
5014 {
5015 // Assign coordinates of dummy element
5016 for (unsigned j = 0; j < 8; j++)
5017 {
5021 switch (j)
5022 {
5023 case 0:
5025 nod_pt->set_value(0, s_tet[0]);
5026 nod_pt->set_value(1, s_tet[1]);
5027 nod_pt->set_value(2, s_tet[2]);
5029 nod_pt->x(0) = x_tet[0];
5030 nod_pt->x(1) = x_tet[1];
5031 nod_pt->x(2) = x_tet[2];
5032 break;
5033 case 1:
5035 nod_pt->set_value(0, s_tet[0]);
5036 nod_pt->set_value(1, s_tet[1]);
5037 nod_pt->set_value(2, s_tet[2]);
5039 nod_pt->x(0) = x_tet[0];
5040 nod_pt->x(1) = x_tet[1];
5041 nod_pt->x(2) = x_tet[2];
5042 break;
5043 case 2:
5045 nod_pt->set_value(0, s_tet[0]);
5046 nod_pt->set_value(1, s_tet[1]);
5047 nod_pt->set_value(2, s_tet[2]);
5049 nod_pt->x(0) = x_tet[0];
5050 nod_pt->x(1) = x_tet[1];
5051 nod_pt->x(2) = x_tet[2];
5052 break;
5053 case 3:
5054 // label 13 in initial sketch: Mid face node on face spanned by
5055 // tet nodes 0,1,3
5056 s_tet[0] = 1.0 / 3.0;
5057 s_tet[1] = 1.0 / 3.0;
5058 s_tet[2] = 0.0;
5059 nod_pt->set_value(0, s_tet[0]);
5060 nod_pt->set_value(1, s_tet[1]);
5061 nod_pt->set_value(2, s_tet[2]);
5063 nod_pt->x(0) = x_tet[0];
5064 nod_pt->x(1) = x_tet[1];
5065 nod_pt->x(2) = x_tet[2];
5066 break;
5067 case 4:
5069 nod_pt->set_value(0, s_tet[0]);
5070 nod_pt->set_value(1, s_tet[1]);
5071 nod_pt->set_value(2, s_tet[2]);
5073 nod_pt->x(0) = x_tet[0];
5074 nod_pt->x(1) = x_tet[1];
5075 nod_pt->x(2) = x_tet[2];
5076 break;
5077 case 5:
5078 // label 11 in initial sketch: Mid face node on face spanned
5079 // by tet nodes 0,1,2
5080 s_tet[0] = 1.0 / 3.0;
5081 s_tet[1] = 1.0 / 3.0;
5082 s_tet[2] = 1.0 / 3.0;
5083 nod_pt->set_value(0, s_tet[0]);
5084 nod_pt->set_value(1, s_tet[1]);
5085 nod_pt->set_value(2, s_tet[2]);
5087 nod_pt->x(0) = x_tet[0];
5088 nod_pt->x(1) = x_tet[1];
5089 nod_pt->x(2) = x_tet[2];
5090 break;
5091 case 6:
5092 // label 12 in initial sketch: Mid face node on face
5093 // spanned by tet nodes 0,2,3
5094 s_tet[0] = 1.0 / 3.0;
5095 s_tet[1] = 0.0;
5096 s_tet[2] = 1.0 / 3.0;
5097 nod_pt->set_value(0, s_tet[0]);
5098 nod_pt->set_value(1, s_tet[1]);
5099 nod_pt->set_value(2, s_tet[2]);
5101 nod_pt->x(0) = x_tet[0];
5102 nod_pt->x(1) = x_tet[1];
5103 nod_pt->x(2) = x_tet[2];
5104 break;
5105 case 7:
5106 // label 14 in initial sketch: Centroid
5107 s_tet[0] = 0.25;
5108 s_tet[1] = 0.25;
5109 s_tet[2] = 0.25;
5110 nod_pt->set_value(0, s_tet[0]);
5111 nod_pt->set_value(1, s_tet[1]);
5112 nod_pt->set_value(2, s_tet[2]);
5114 nod_pt->x(0) = x_tet[0];
5115 nod_pt->x(1) = x_tet[1];
5116 nod_pt->x(2) = x_tet[2];
5117 break;
5118 }
5119 }
5120
5121 // Create actual zeroth brick element
5122 FiniteElement* el_pt = new ELEMENT;
5124 Element_pt.push_back(el_pt);
5125
5126 TFace face0(
5128
5129 TFace face1(
5131
5132 TFace face2(
5134
5135 // Tet vertex nodes along edges emanating from node 0 in brick
5137 tet_edge_node[0].resize(2);
5138 tet_edge_node[0][0] = 4;
5139 tet_edge_node[0][1] = 1;
5140 tet_edge_node[1].resize(2);
5141 tet_edge_node[1][0] = 6;
5142 tet_edge_node[1][1] = 3;
5143 tet_edge_node[2].resize(2);
5144 tet_edge_node[2][0] = 5;
5145 tet_edge_node[2][1] = 2;
5146
5147 // Node number of tet vertex that node 0 in brick is centred on
5148 unsigned central_tet_vertex = 0;
5149
5150 Node* tet_node_pt = 0;
5151 Node* old_node_pt = 0;
5152
5153 // Corner node
5154 {
5155 unsigned j = 0;
5156
5157 // Need new node?
5160 if (old_node_pt == 0)
5161 {
5162 Node* new_node_pt = 0;
5163 if (tet_node_pt->is_on_boundary())
5164 {
5166 }
5167 else
5168 {
5170 }
5172 Node_pt.push_back(new_node_pt);
5173 Vector<double> s(3);
5177 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5179 new_node_pt->x(0) = x_tet[0];
5180 new_node_pt->x(1) = x_tet[1];
5181 new_node_pt->x(2) = x_tet[2];
5182 }
5183 // Node already exists
5184 else
5185 {
5187 }
5188 }
5189
5190 // Brick vertex node coindides with mid-edge node on tet edge 0
5191 {
5192 unsigned j = 2;
5193
5194 // Need new node?
5197 if (old_node_pt == 0)
5198 {
5199 Node* new_node_pt = 0;
5200 if (tet_node_pt->is_on_boundary())
5201 {
5203 }
5204 else
5205 {
5207 }
5209 Node_pt.push_back(new_node_pt);
5210 Vector<double> s(3);
5214 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5216 new_node_pt->x(0) = x_tet[0];
5217 new_node_pt->x(1) = x_tet[1];
5218 new_node_pt->x(2) = x_tet[2];
5219 }
5220 // Node already exists
5221 else
5222 {
5224 }
5225 }
5226
5227 // Brick vertex node coindides with mid vertex node of tet edge 1
5228 {
5229 unsigned j = 6;
5230
5231 // Need new node?
5234 if (old_node_pt == 0)
5235 {
5236 Node* new_node_pt = 0;
5237 if (tet_node_pt->is_on_boundary())
5238 {
5240 }
5241 else
5242 {
5244 }
5246 Node_pt.push_back(new_node_pt);
5247 Vector<double> s(3);
5251 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5253 new_node_pt->x(0) = x_tet[0];
5254 new_node_pt->x(1) = x_tet[1];
5255 new_node_pt->x(2) = x_tet[2];
5256 }
5257 // Node already exists
5258 else
5259 {
5261 }
5262 }
5263
5264 // Brick vertex node coindides with mid-vertex node of tet edge 2
5265 {
5266 unsigned j = 18;
5267
5268 // Need new node?
5271 if (old_node_pt == 0)
5272 {
5273 Node* new_node_pt = 0;
5274 if (tet_node_pt->is_on_boundary())
5275 {
5277 }
5278 else
5279 {
5281 }
5283 Node_pt.push_back(new_node_pt);
5284 Vector<double> s(3);
5288 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5290 new_node_pt->x(0) = x_tet[0];
5291 new_node_pt->x(1) = x_tet[1];
5292 new_node_pt->x(2) = x_tet[2];
5293 }
5294 // Node already exists
5295 else
5296 {
5298 }
5299 }
5300
5301 // Brick vertex node in the middle of tet face0, spanned by
5302 // tet vertices 0, 1, 2. Enumerated "11" in initial sketch.
5303 {
5304 unsigned j = 20;
5305
5306 // Need new node?
5308 if (old_node_pt == 0)
5309 {
5310 Node* new_node_pt = 0;
5311 if (face0.is_boundary_face())
5312 {
5314 }
5315 else
5316 {
5318 }
5320 Node_pt.push_back(new_node_pt);
5321 Vector<double> s(3);
5325 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5327 new_node_pt->x(0) = x_tet[0];
5328 new_node_pt->x(1) = x_tet[1];
5329 new_node_pt->x(2) = x_tet[2];
5330 }
5331 // Node already exists
5332 else
5333 {
5335 }
5336 }
5337
5338 // Brick vertex node in the middle of tet face1, spanned by
5339 // tet vertices 0, 2, 3. Enumerated "12" in initial sketch.
5340 {
5341 unsigned j = 24;
5342
5343 // Need new node?
5345 if (old_node_pt == 0)
5346 {
5347 Node* new_node_pt = 0;
5348 if (face1.is_boundary_face())
5349 {
5351 }
5352 else
5353 {
5355 }
5357 Node_pt.push_back(new_node_pt);
5358 Vector<double> s(3);
5362 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5364 new_node_pt->x(0) = x_tet[0];
5365 new_node_pt->x(1) = x_tet[1];
5366 new_node_pt->x(2) = x_tet[2];
5367 }
5368 // Node already exists
5369 else
5370 {
5372 }
5373 }
5374
5375 // Brick vertex node in the middle of tet face2, spanned by
5376 // tet vertices 0, 1, 3. Enumerated "13" in initial sketch.
5377 {
5378 unsigned j = 8;
5379
5380 // Need new node?
5382 if (old_node_pt == 0)
5383 {
5384 Node* new_node_pt = 0;
5385 if (face2.is_boundary_face())
5386 {
5388 }
5389 else
5390 {
5392 }
5394 Node_pt.push_back(new_node_pt);
5395 Vector<double> s(3);
5399 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5401 new_node_pt->x(0) = x_tet[0];
5402 new_node_pt->x(1) = x_tet[1];
5403 new_node_pt->x(2) = x_tet[2];
5404 }
5405 // Node already exists
5406 else
5407 {
5409 }
5410 }
5411
5412 // Brick vertex node in centroid of tet. Only built for first element.
5413 // Enumerated "13" in initial sketch.
5414 {
5415 unsigned j = 26;
5416
5417 // Always new
5418 {
5421 Node_pt.push_back(new_node_pt);
5422 Vector<double> s(3);
5426 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5428 new_node_pt->x(0) = x_tet[0];
5429 new_node_pt->x(1) = x_tet[1];
5430 new_node_pt->x(2) = x_tet[2];
5431 }
5432 }
5433
5434 // Internal brick node -- always built
5435 {
5436 unsigned j = 13;
5437
5438 // Always new
5439 {
5441 Node_pt.push_back(new_node_pt);
5442 Vector<double> s(3);
5446 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5448 new_node_pt->x(0) = x_tet[0];
5449 new_node_pt->x(1) = x_tet[1];
5450 new_node_pt->x(2) = x_tet[2];
5451 }
5452 }
5453
5454 // Brick edge node between brick nodes 0 and 2
5455 {
5456 unsigned j = 1;
5457
5458 // Need new node?
5459 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
5461 if (old_node_pt == 0)
5462 {
5463 Node* new_node_pt = 0;
5464 if (edge.is_boundary_edge())
5465 {
5467 }
5468 else
5469 {
5471 }
5473 Node_pt.push_back(new_node_pt);
5474 Vector<double> s(3);
5478 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5480 new_node_pt->x(0) = x_tet[0];
5481 new_node_pt->x(1) = x_tet[1];
5482 new_node_pt->x(2) = x_tet[2];
5483 }
5484 // Node already exists
5485 else
5486 {
5488 }
5489 }
5490
5491 // Brick edge node between brick nodes 0 and 6
5492 {
5493 unsigned j = 3;
5494
5495 // Need new node?
5496 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
5498 if (old_node_pt == 0)
5499 {
5500 Node* new_node_pt = 0;
5501 if (edge.is_boundary_edge())
5502 {
5504 }
5505 else
5506 {
5508 }
5510 Node_pt.push_back(new_node_pt);
5511 Vector<double> s(3);
5515 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5517 new_node_pt->x(0) = x_tet[0];
5518 new_node_pt->x(1) = x_tet[1];
5519 new_node_pt->x(2) = x_tet[2];
5520 }
5521 // Node already exists
5522 else
5523 {
5525 }
5526 }
5527
5528 // Brick edge node between brick nodes 2 and 8
5529 {
5530 unsigned j = 5;
5531
5532 // Need new node?
5533 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
5535 if (old_node_pt == 0)
5536 {
5537 Node* new_node_pt = 0;
5538 if (edge.is_boundary_edge())
5539 {
5541 }
5542 else
5543 {
5545 }
5547 Node_pt.push_back(new_node_pt);
5548 Vector<double> s(3);
5552 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5554 new_node_pt->x(0) = x_tet[0];
5555 new_node_pt->x(1) = x_tet[1];
5556 new_node_pt->x(2) = x_tet[2];
5557 }
5558 // Node already exists
5559 else
5560 {
5562 }
5563 }
5564
5565 // Brick edge node between brick nodes 6 and 8
5566 {
5567 unsigned j = 7;
5568
5569 // Need new node?
5570 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
5572 if (old_node_pt == 0)
5573 {
5574 Node* new_node_pt = 0;
5575 if (edge.is_boundary_edge())
5576 {
5578 }
5579 else
5580 {
5582 }
5584 Node_pt.push_back(new_node_pt);
5585 Vector<double> s(3);
5589 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5591 new_node_pt->x(0) = x_tet[0];
5592 new_node_pt->x(1) = x_tet[1];
5593 new_node_pt->x(2) = x_tet[2];
5594 }
5595 // Node already exists
5596 else
5597 {
5599 }
5600 }
5601
5602 // Brick edge node between brick nodes 18 and 20
5603 {
5604 unsigned j = 19;
5605
5606 // Need new node?
5607 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
5609 if (old_node_pt == 0)
5610 {
5611 Node* new_node_pt = 0;
5612 if (edge.is_boundary_edge())
5613 {
5615 }
5616 else
5617 {
5619 }
5621 Node_pt.push_back(new_node_pt);
5622 Vector<double> s(3);
5626 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5628 new_node_pt->x(0) = x_tet[0];
5629 new_node_pt->x(1) = x_tet[1];
5630 new_node_pt->x(2) = x_tet[2];
5631 }
5632 // Node already exists
5633 else
5634 {
5636 }
5637 }
5638
5639 // Brick edge node between brick nodes 18 and 24
5640 {
5641 unsigned j = 21;
5642
5643 // Need new node?
5644 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
5646 if (old_node_pt == 0)
5647 {
5648 Node* new_node_pt = 0;
5649 if (edge.is_boundary_edge())
5650 {
5652 }
5653 else
5654 {
5656 }
5658 Node_pt.push_back(new_node_pt);
5659 Vector<double> s(3);
5663 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5665 new_node_pt->x(0) = x_tet[0];
5666 new_node_pt->x(1) = x_tet[1];
5667 new_node_pt->x(2) = x_tet[2];
5668 }
5669 // Node already exists
5670 else
5671 {
5673 }
5674 }
5675
5676 // Brick edge node between brick nodes 20 and 26
5677 {
5678 unsigned j = 23;
5679
5680 // Need new node?
5681 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
5683 if (old_node_pt == 0)
5684 {
5685 Node* new_node_pt = 0;
5686 if (edge.is_boundary_edge())
5687 {
5689 }
5690 else
5691 {
5693 }
5695 Node_pt.push_back(new_node_pt);
5696 Vector<double> s(3);
5700 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5702 new_node_pt->x(0) = x_tet[0];
5703 new_node_pt->x(1) = x_tet[1];
5704 new_node_pt->x(2) = x_tet[2];
5705 }
5706 // Node already exists
5707 else
5708 {
5710 }
5711 }
5712
5713 // Brick edge node between brick nodes 24 and 26
5714 {
5715 unsigned j = 25;
5716
5717 // Need new node?
5718 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
5720 if (old_node_pt == 0)
5721 {
5722 Node* new_node_pt = 0;
5723 if (edge.is_boundary_edge())
5724 {
5726 }
5727 else
5728 {
5730 }
5732 Node_pt.push_back(new_node_pt);
5733 Vector<double> s(3);
5737 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5739 new_node_pt->x(0) = x_tet[0];
5740 new_node_pt->x(1) = x_tet[1];
5741 new_node_pt->x(2) = x_tet[2];
5742 }
5743 // Node already exists
5744 else
5745 {
5747 }
5748 }
5749
5750 // Brick edge node between brick nodes 0 and 18
5751 {
5752 unsigned j = 9;
5753
5754 // Need new node?
5755 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
5757 if (old_node_pt == 0)
5758 {
5759 Node* new_node_pt = 0;
5760 if (edge.is_boundary_edge())
5761 {
5763 }
5764 else
5765 {
5767 }
5769 Node_pt.push_back(new_node_pt);
5770 Vector<double> s(3);
5774 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5776 new_node_pt->x(0) = x_tet[0];
5777 new_node_pt->x(1) = x_tet[1];
5778 new_node_pt->x(2) = x_tet[2];
5779 }
5780 // Node already exists
5781 else
5782 {
5784 }
5785 }
5786
5787 // Brick edge node between brick nodes 2 and 20
5788 {
5789 unsigned j = 11;
5790
5791 // Need new node?
5792 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
5794 if (old_node_pt == 0)
5795 {
5796 Node* new_node_pt = 0;
5797 if (edge.is_boundary_edge())
5798 {
5800 }
5801 else
5802 {
5804 }
5806 Node_pt.push_back(new_node_pt);
5807 Vector<double> s(3);
5811 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5813 new_node_pt->x(0) = x_tet[0];
5814 new_node_pt->x(1) = x_tet[1];
5815 new_node_pt->x(2) = x_tet[2];
5816 }
5817 // Node already exists
5818 else
5819 {
5821 }
5822 }
5823
5824 // Brick edge node between brick nodes 6 and 24
5825 {
5826 unsigned j = 15;
5827
5828 // Need new node?
5829 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
5831 if (old_node_pt == 0)
5832 {
5833 Node* new_node_pt = 0;
5834 if (edge.is_boundary_edge())
5835 {
5837 }
5838 else
5839 {
5841 }
5843 Node_pt.push_back(new_node_pt);
5844 Vector<double> s(3);
5848 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5850 new_node_pt->x(0) = x_tet[0];
5851 new_node_pt->x(1) = x_tet[1];
5852 new_node_pt->x(2) = x_tet[2];
5853 }
5854 // Node already exists
5855 else
5856 {
5858 }
5859 }
5860
5861 // Brick edge node between brick nodes 8 and 26
5862 {
5863 unsigned j = 17;
5864
5865 // Need new node?
5866 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
5868 if (old_node_pt == 0)
5869 {
5870 Node* new_node_pt = 0;
5871 if (edge.is_boundary_edge())
5872 {
5874 }
5875 else
5876 {
5878 }
5880 Node_pt.push_back(new_node_pt);
5881 Vector<double> s(3);
5885 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5887 new_node_pt->x(0) = x_tet[0];
5888 new_node_pt->x(1) = x_tet[1];
5889 new_node_pt->x(2) = x_tet[2];
5890 }
5891 // Node already exists
5892 else
5893 {
5895 }
5896 }
5897
5898 // Mid brick-face node associated with face
5899 // spanned by mid-vertex nodes associated with tet edges 0 and 2
5900 {
5901 unsigned j = 10;
5902
5903 // Need new node?
5907
5909 if (old_node_pt == 0)
5910 {
5911 Node* new_node_pt = 0;
5912 if (face.is_boundary_face())
5913 {
5915 }
5916 else
5917 {
5919 }
5921 Node_pt.push_back(new_node_pt);
5922 Vector<double> s(3);
5926 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5928 new_node_pt->x(0) = x_tet[0];
5929 new_node_pt->x(1) = x_tet[1];
5930 new_node_pt->x(2) = x_tet[2];
5931 }
5932 // Node already exists
5933 else
5934 {
5936 }
5937 }
5938
5939 // Mid brick-face node associated with face
5940 // spanned by mid-vertex nodes associated with tet edges 1 and 2
5941 {
5942 unsigned j = 12;
5943
5944 // Need new node?
5948
5950 if (old_node_pt == 0)
5951 {
5952 Node* new_node_pt = 0;
5953 if (face.is_boundary_face())
5954 {
5956 }
5957 else
5958 {
5960 }
5962 Node_pt.push_back(new_node_pt);
5963 Vector<double> s(3);
5967 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
5969 new_node_pt->x(0) = x_tet[0];
5970 new_node_pt->x(1) = x_tet[1];
5971 new_node_pt->x(2) = x_tet[2];
5972 }
5973 // Node already exists
5974 else
5975 {
5977 }
5978 }
5979
5980 // Mid brick-face node associated with face
5981 // spanned by mid-vertex nodes associated with tet edges 0 and 1
5982 {
5983 unsigned j = 4;
5984
5985 // Need new node?
5989
5991 if (old_node_pt == 0)
5992 {
5993 Node* new_node_pt = 0;
5994 if (face.is_boundary_face())
5995 {
5997 }
5998 else
5999 {
6001 }
6003 Node_pt.push_back(new_node_pt);
6004 Vector<double> s(3);
6008 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
6010 new_node_pt->x(0) = x_tet[0];
6011 new_node_pt->x(1) = x_tet[1];
6012 new_node_pt->x(2) = x_tet[2];
6013 }
6014 // Node already exists
6015 else
6016 {
6018 }
6019 }
6020
6021 // Top mid brick-face node -- only built by first element
6022 {
6023 unsigned j = 22;
6025 Node_pt.push_back(new_node_pt);
6026 Vector<double> s(3);
6030 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
6033 new_node_pt->x(0) = x_tet[0];
6034 new_node_pt->x(1) = x_tet[1];
6035 new_node_pt->x(2) = x_tet[2];
6036 }
6037
6038 // Right mid brick-face node -- only built by first element
6039 {
6040 unsigned j = 14;
6042 Node_pt.push_back(new_node_pt);
6043 Vector<double> s(3);
6047 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
6050 new_node_pt->x(0) = x_tet[0];
6051 new_node_pt->x(1) = x_tet[1];
6052 new_node_pt->x(2) = x_tet[2];
6053 }
6054
6055 // Back mid brick-face node -- only built by first element
6056 {
6057 unsigned j = 16;
6059 Node_pt.push_back(new_node_pt);
6060 Vector<double> s(3);
6064 dummy_q_el_pt[0]->interpolated_s_tet(s, s_tet);
6067 new_node_pt->x(0) = x_tet[0];
6068 new_node_pt->x(1) = x_tet[1];
6069 new_node_pt->x(2) = x_tet[2];
6070 }
6071 }
6072
6073 // Second brick element is centred at node 1 of tet:
6074 //--------------------------------------------------
6075 {
6076 // Assign coordinates of dummy element
6077 for (unsigned j = 0; j < 8; j++)
6078 {
6082 switch (j)
6083 {
6084 case 0:
6086 nod_pt->set_value(0, s_tet[0]);
6087 nod_pt->set_value(1, s_tet[1]);
6088 nod_pt->set_value(2, s_tet[2]);
6090 nod_pt->x(0) = x_tet[0];
6091 nod_pt->x(1) = x_tet[1];
6092 nod_pt->x(2) = x_tet[2];
6093 break;
6094 case 1:
6096 nod_pt->set_value(0, s_tet[0]);
6097 nod_pt->set_value(1, s_tet[1]);
6098 nod_pt->set_value(2, s_tet[2]);
6100 nod_pt->x(0) = x_tet[0];
6101 nod_pt->x(1) = x_tet[1];
6102 nod_pt->x(2) = x_tet[2];
6103 break;
6104 case 2:
6106 nod_pt->set_value(0, s_tet[0]);
6107 nod_pt->set_value(1, s_tet[1]);
6108 nod_pt->set_value(2, s_tet[2]);
6110 nod_pt->x(0) = x_tet[0];
6111 nod_pt->x(1) = x_tet[1];
6112 nod_pt->x(2) = x_tet[2];
6113 break;
6114 case 3:
6115 // label 13 in initial sketch: Mid face node on face
6116 // spanned by tet nodes 0,1,3
6117 s_tet[0] = 1.0 / 3.0;
6118 s_tet[1] = 1.0 / 3.0;
6119 s_tet[2] = 0.0;
6120 nod_pt->set_value(0, s_tet[0]);
6121 nod_pt->set_value(1, s_tet[1]);
6122 nod_pt->set_value(2, s_tet[2]);
6124 nod_pt->x(0) = x_tet[0];
6125 nod_pt->x(1) = x_tet[1];
6126 nod_pt->x(2) = x_tet[2];
6127 break;
6128 case 4:
6130 nod_pt->set_value(0, s_tet[0]);
6131 nod_pt->set_value(1, s_tet[1]);
6132 nod_pt->set_value(2, s_tet[2]);
6134 nod_pt->x(0) = x_tet[0];
6135 nod_pt->x(1) = x_tet[1];
6136 nod_pt->x(2) = x_tet[2];
6137 break;
6138 case 5:
6139 // label 10 in initial sketch: Mid face node on face
6140 // spanned by tet nodes 1,2,3
6141 s_tet[0] = 0.0;
6142 s_tet[1] = 1.0 / 3.0;
6143 s_tet[2] = 1.0 / 3.0;
6144 nod_pt->set_value(0, s_tet[0]);
6145 nod_pt->set_value(1, s_tet[1]);
6146 nod_pt->set_value(2, s_tet[2]);
6148 nod_pt->x(0) = x_tet[0];
6149 nod_pt->x(1) = x_tet[1];
6150 nod_pt->x(2) = x_tet[2];
6151 break;
6152 case 6:
6153 // label 11 in initial sketch: Mid face node on face
6154 // spanned by tet nodes 0,1,2
6155 s_tet[0] = 1.0 / 3.0;
6156 s_tet[1] = 1.0 / 3.0;
6157 s_tet[2] = 1.0 / 3.0;
6158 nod_pt->set_value(0, s_tet[0]);
6159 nod_pt->set_value(1, s_tet[1]);
6160 nod_pt->set_value(2, s_tet[2]);
6162 nod_pt->x(0) = x_tet[0];
6163 nod_pt->x(1) = x_tet[1];
6164 nod_pt->x(2) = x_tet[2];
6165 break;
6166 case 7:
6167 // label 14 in initial sketch: Centroid
6168 s_tet[0] = 0.25;
6169 s_tet[1] = 0.25;
6170 s_tet[2] = 0.25;
6171 nod_pt->set_value(0, s_tet[0]);
6172 nod_pt->set_value(1, s_tet[1]);
6173 nod_pt->set_value(2, s_tet[2]);
6175 nod_pt->x(0) = x_tet[0];
6176 nod_pt->x(1) = x_tet[1];
6177 nod_pt->x(2) = x_tet[2];
6178 break;
6179 }
6180 }
6181
6182 // Create actual first brick element
6183 FiniteElement* el_pt = new ELEMENT;
6185 Element_pt.push_back(el_pt);
6186
6187 TFace face0(
6189
6190 TFace face1(
6192
6193 TFace face2(
6195
6196 // Tet vertex nodes along edges emanating from node 0 in brick
6198 tet_edge_node[0].resize(2);
6199 tet_edge_node[0][0] = 9;
6200 tet_edge_node[0][1] = 3;
6201 tet_edge_node[1].resize(2);
6202 tet_edge_node[1][0] = 4;
6203 tet_edge_node[1][1] = 0;
6204 tet_edge_node[2].resize(2);
6205 tet_edge_node[2][0] = 7;
6206 tet_edge_node[2][1] = 2;
6207
6208 // Node number of tet vertex that node 0 in brick is centred on
6209 unsigned central_tet_vertex = 1;
6210
6211 Node* tet_node_pt = 0;
6212 Node* old_node_pt = 0;
6213
6214 // Corner node
6215 {
6216 unsigned j = 0;
6217
6218 // Need new node?
6221 if (old_node_pt == 0)
6222 {
6223 Node* new_node_pt = 0;
6224 if (tet_node_pt->is_on_boundary())
6225 {
6227 }
6228 else
6229 {
6231 }
6233 Node_pt.push_back(new_node_pt);
6234 Vector<double> s(3);
6238 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6240 new_node_pt->x(0) = x_tet[0];
6241 new_node_pt->x(1) = x_tet[1];
6242 new_node_pt->x(2) = x_tet[2];
6243 }
6244 // Node already exists
6245 else
6246 {
6248 }
6249 }
6250
6251 // Brick vertex node coindides with mid-edge node on tet edge 0
6252 {
6253 unsigned j = 2;
6254
6255 // Need new node?
6258 if (old_node_pt == 0)
6259 {
6260 Node* new_node_pt = 0;
6261 if (tet_node_pt->is_on_boundary())
6262 {
6264 }
6265 else
6266 {
6268 }
6270 Node_pt.push_back(new_node_pt);
6271 Vector<double> s(3);
6275 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6277 new_node_pt->x(0) = x_tet[0];
6278 new_node_pt->x(1) = x_tet[1];
6279 new_node_pt->x(2) = x_tet[2];
6280 }
6281 // Node already exists
6282 else
6283 {
6285 }
6286 }
6287
6288 // Brick vertex node coindides with mid vertex node of tet edge 1
6289 {
6290 unsigned j = 6;
6291
6292 // Need new node?
6295 if (old_node_pt == 0)
6296 {
6297 Node* new_node_pt = 0;
6298 if (tet_node_pt->is_on_boundary())
6299 {
6301 }
6302 else
6303 {
6305 }
6307 Node_pt.push_back(new_node_pt);
6308 Vector<double> s(3);
6312 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6314 new_node_pt->x(0) = x_tet[0];
6315 new_node_pt->x(1) = x_tet[1];
6316 new_node_pt->x(2) = x_tet[2];
6317 }
6318 // Node already exists
6319 else
6320 {
6322 }
6323 }
6324
6325 // Brick vertex node coindides with mid-vertex node of tet edge 2
6326 {
6327 unsigned j = 18;
6328
6329 // Need new node?
6332 if (old_node_pt == 0)
6333 {
6334 Node* new_node_pt = 0;
6335 if (tet_node_pt->is_on_boundary())
6336 {
6338 }
6339 else
6340 {
6342 }
6344 Node_pt.push_back(new_node_pt);
6345 Vector<double> s(3);
6349 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6351 new_node_pt->x(0) = x_tet[0];
6352 new_node_pt->x(1) = x_tet[1];
6353 new_node_pt->x(2) = x_tet[2];
6354 }
6355 // Node already exists
6356 else
6357 {
6359 }
6360 }
6361
6362 // Brick vertex node in the middle of tet face0
6363 {
6364 unsigned j = 20;
6365
6366 // Need new node?
6368 if (old_node_pt == 0)
6369 {
6370 Node* new_node_pt = 0;
6371 if (face0.is_boundary_face())
6372 {
6374 }
6375 else
6376 {
6378 }
6380 Node_pt.push_back(new_node_pt);
6381 Vector<double> s(3);
6385 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6387 new_node_pt->x(0) = x_tet[0];
6388 new_node_pt->x(1) = x_tet[1];
6389 new_node_pt->x(2) = x_tet[2];
6390 }
6391 // Node already exists
6392 else
6393 {
6395 }
6396 }
6397
6398 // Brick vertex node in the middle of tet face1
6399 {
6400 unsigned j = 24;
6401
6402 // Need new node?
6404 if (old_node_pt == 0)
6405 {
6406 Node* new_node_pt = 0;
6407 if (face1.is_boundary_face())
6408 {
6410 }
6411 else
6412 {
6414 }
6416 Node_pt.push_back(new_node_pt);
6417 Vector<double> s(3);
6421 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6423 new_node_pt->x(0) = x_tet[0];
6424 new_node_pt->x(1) = x_tet[1];
6425 new_node_pt->x(2) = x_tet[2];
6426 }
6427 // Node already exists
6428 else
6429 {
6431 }
6432 }
6433
6434 // Brick vertex node in the middle of face2
6435 {
6436 unsigned j = 8;
6437
6438 // Need new node?
6440 if (old_node_pt == 0)
6441 {
6442 Node* new_node_pt = 0;
6443 if (face2.is_boundary_face())
6444 {
6446 }
6447 else
6448 {
6450 }
6452 Node_pt.push_back(new_node_pt);
6453 Vector<double> s(3);
6457 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6459 new_node_pt->x(0) = x_tet[0];
6460 new_node_pt->x(1) = x_tet[1];
6461 new_node_pt->x(2) = x_tet[2];
6462 }
6463 // Node already exists
6464 else
6465 {
6467 }
6468 }
6469
6470 // Brick vertex node in centroid of tet. Only built for first element.
6471 // Enumerated "13" in initial sketch.
6472 {
6473 unsigned j = 26;
6474
6475 // Always copied
6477 }
6478
6479 // Internal brick node -- always built
6480 {
6481 unsigned j = 13;
6482
6483 // Always new
6484 {
6486 Node_pt.push_back(new_node_pt);
6487 Vector<double> s(3);
6491 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6493 new_node_pt->x(0) = x_tet[0];
6494 new_node_pt->x(1) = x_tet[1];
6495 new_node_pt->x(2) = x_tet[2];
6496 }
6497 }
6498
6499 // Brick edge node between brick nodes 0 and 2
6500 {
6501 unsigned j = 1;
6502
6503 // Need new node?
6504 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
6506 if (old_node_pt == 0)
6507 {
6508 Node* new_node_pt = 0;
6509 if (edge.is_boundary_edge())
6510 {
6512 }
6513 else
6514 {
6516 }
6518 Node_pt.push_back(new_node_pt);
6519 Vector<double> s(3);
6523 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6525 new_node_pt->x(0) = x_tet[0];
6526 new_node_pt->x(1) = x_tet[1];
6527 new_node_pt->x(2) = x_tet[2];
6528 }
6529 // Node already exists
6530 else
6531 {
6533 }
6534 }
6535
6536 // Brick edge node between brick nodes 0 and 6
6537 {
6538 unsigned j = 3;
6539
6540 // Need new node?
6541 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
6543 if (old_node_pt == 0)
6544 {
6545 Node* new_node_pt = 0;
6546 if (edge.is_boundary_edge())
6547 {
6549 }
6550 else
6551 {
6553 }
6555 Node_pt.push_back(new_node_pt);
6556 Vector<double> s(3);
6560 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6562 new_node_pt->x(0) = x_tet[0];
6563 new_node_pt->x(1) = x_tet[1];
6564 new_node_pt->x(2) = x_tet[2];
6565 }
6566 // Node already exists
6567 else
6568 {
6570 }
6571 }
6572
6573 // Brick edge node between brick nodes 2 and 8
6574 {
6575 unsigned j = 5;
6576
6577 // Need new node?
6578 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
6580 if (old_node_pt == 0)
6581 {
6582 Node* new_node_pt = 0;
6583 if (edge.is_boundary_edge())
6584 {
6586 }
6587 else
6588 {
6590 }
6592 Node_pt.push_back(new_node_pt);
6593 Vector<double> s(3);
6597 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6599 new_node_pt->x(0) = x_tet[0];
6600 new_node_pt->x(1) = x_tet[1];
6601 new_node_pt->x(2) = x_tet[2];
6602 }
6603 // Node already exists
6604 else
6605 {
6607 }
6608 }
6609
6610 // Brick edge node between brick nodes 6 and 8
6611 {
6612 unsigned j = 7;
6613
6614 // Need new node?
6615 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
6617 if (old_node_pt == 0)
6618 {
6619 Node* new_node_pt = 0;
6620 if (edge.is_boundary_edge())
6621 {
6623 }
6624 else
6625 {
6627 }
6629 Node_pt.push_back(new_node_pt);
6630 Vector<double> s(3);
6634 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6636 new_node_pt->x(0) = x_tet[0];
6637 new_node_pt->x(1) = x_tet[1];
6638 new_node_pt->x(2) = x_tet[2];
6639 }
6640 // Node already exists
6641 else
6642 {
6644 }
6645 }
6646
6647 // Brick edge node between brick nodes 18 and 20
6648 {
6649 unsigned j = 19;
6650
6651 // Need new node?
6652 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
6654 if (old_node_pt == 0)
6655 {
6656 Node* new_node_pt = 0;
6657 if (edge.is_boundary_edge())
6658 {
6660 }
6661 else
6662 {
6664 }
6666 Node_pt.push_back(new_node_pt);
6667 Vector<double> s(3);
6671 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6673 new_node_pt->x(0) = x_tet[0];
6674 new_node_pt->x(1) = x_tet[1];
6675 new_node_pt->x(2) = x_tet[2];
6676 }
6677 // Node already exists
6678 else
6679 {
6681 }
6682 }
6683
6684 // Brick edge node between brick nodes 18 and 24
6685 {
6686 unsigned j = 21;
6687
6688 // Need new node?
6689 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
6691 if (old_node_pt == 0)
6692 {
6693 Node* new_node_pt = 0;
6694 if (edge.is_boundary_edge())
6695 {
6697 }
6698 else
6699 {
6701 }
6703 Node_pt.push_back(new_node_pt);
6704 Vector<double> s(3);
6708 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6710 new_node_pt->x(0) = x_tet[0];
6711 new_node_pt->x(1) = x_tet[1];
6712 new_node_pt->x(2) = x_tet[2];
6713 }
6714 // Node already exists
6715 else
6716 {
6718 }
6719 }
6720
6721 // Brick edge node between brick nodes 20 and 26
6722 {
6723 unsigned j = 23;
6724
6725 // Need new node?
6726 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
6728 if (old_node_pt == 0)
6729 {
6730 Node* new_node_pt = 0;
6731 if (edge.is_boundary_edge())
6732 {
6734 }
6735 else
6736 {
6738 }
6740 Node_pt.push_back(new_node_pt);
6741 Vector<double> s(3);
6745 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6747 new_node_pt->x(0) = x_tet[0];
6748 new_node_pt->x(1) = x_tet[1];
6749 new_node_pt->x(2) = x_tet[2];
6750 }
6751 // Node already exists
6752 else
6753 {
6755 }
6756 }
6757
6758 // Brick edge node between brick nodes 24 and 26
6759 {
6760 unsigned j = 25;
6761
6762 // Need new node?
6763 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
6765 if (old_node_pt == 0)
6766 {
6767 Node* new_node_pt = 0;
6768 if (edge.is_boundary_edge())
6769 {
6771 }
6772 else
6773 {
6775 }
6777 Node_pt.push_back(new_node_pt);
6778 Vector<double> s(3);
6782 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6784 new_node_pt->x(0) = x_tet[0];
6785 new_node_pt->x(1) = x_tet[1];
6786 new_node_pt->x(2) = x_tet[2];
6787 }
6788 // Node already exists
6789 else
6790 {
6792 }
6793 }
6794
6795 // Brick edge node between brick nodes 0 and 18
6796 {
6797 unsigned j = 9;
6798
6799 // Need new node?
6800 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
6802 if (old_node_pt == 0)
6803 {
6804 Node* new_node_pt = 0;
6805 if (edge.is_boundary_edge())
6806 {
6808 }
6809 else
6810 {
6812 }
6814 Node_pt.push_back(new_node_pt);
6815 Vector<double> s(3);
6819 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6821 new_node_pt->x(0) = x_tet[0];
6822 new_node_pt->x(1) = x_tet[1];
6823 new_node_pt->x(2) = x_tet[2];
6824 }
6825 // Node already exists
6826 else
6827 {
6829 }
6830 }
6831
6832 // Brick edge node between brick nodes 2 and 20
6833 {
6834 unsigned j = 11;
6835
6836 // Need new node?
6837 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
6839 if (old_node_pt == 0)
6840 {
6841 Node* new_node_pt = 0;
6842 if (edge.is_boundary_edge())
6843 {
6845 }
6846 else
6847 {
6849 }
6851 Node_pt.push_back(new_node_pt);
6852 Vector<double> s(3);
6856 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6858 new_node_pt->x(0) = x_tet[0];
6859 new_node_pt->x(1) = x_tet[1];
6860 new_node_pt->x(2) = x_tet[2];
6861 }
6862 // Node already exists
6863 else
6864 {
6866 }
6867 }
6868
6869 // Brick edge node between brick nodes 6 and 24
6870 {
6871 unsigned j = 15;
6872
6873 // Need new node?
6874 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
6876 if (old_node_pt == 0)
6877 {
6878 Node* new_node_pt = 0;
6879 if (edge.is_boundary_edge())
6880 {
6882 }
6883 else
6884 {
6886 }
6888 Node_pt.push_back(new_node_pt);
6889 Vector<double> s(3);
6893 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6895 new_node_pt->x(0) = x_tet[0];
6896 new_node_pt->x(1) = x_tet[1];
6897 new_node_pt->x(2) = x_tet[2];
6898 }
6899 // Node already exists
6900 else
6901 {
6903 }
6904 }
6905
6906 // Brick edge node between brick nodes 8 and 26
6907 {
6908 unsigned j = 17;
6909
6910 // Need new node?
6911 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
6913 if (old_node_pt == 0)
6914 {
6915 Node* new_node_pt = 0;
6916 if (edge.is_boundary_edge())
6917 {
6919 }
6920 else
6921 {
6923 }
6925 Node_pt.push_back(new_node_pt);
6926 Vector<double> s(3);
6930 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6932 new_node_pt->x(0) = x_tet[0];
6933 new_node_pt->x(1) = x_tet[1];
6934 new_node_pt->x(2) = x_tet[2];
6935 }
6936 // Node already exists
6937 else
6938 {
6940 }
6941 }
6942
6943 // Mid brick-face node associated with face
6944 // spanned by mid-vertex nodes associated with tet edges 0 and 2
6945 {
6946 unsigned j = 10;
6947
6948 // Need new node?
6952
6954 if (old_node_pt == 0)
6955 {
6956 Node* new_node_pt = 0;
6957 if (face.is_boundary_face())
6958 {
6960 }
6961 else
6962 {
6964 }
6966 Node_pt.push_back(new_node_pt);
6967 Vector<double> s(3);
6971 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
6973 new_node_pt->x(0) = x_tet[0];
6974 new_node_pt->x(1) = x_tet[1];
6975 new_node_pt->x(2) = x_tet[2];
6976 }
6977 // Node already exists
6978 else
6979 {
6981 }
6982 }
6983
6984 // Mid brick-face node associated with face
6985 // spanned by mid-vertex nodes associated with tet edges 1 and 2
6986 {
6987 unsigned j = 12;
6988
6989 // Need new node?
6993
6995 if (old_node_pt == 0)
6996 {
6997 Node* new_node_pt = 0;
6998 if (face.is_boundary_face())
6999 {
7001 }
7002 else
7003 {
7005 }
7007 Node_pt.push_back(new_node_pt);
7008 Vector<double> s(3);
7012 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
7014 new_node_pt->x(0) = x_tet[0];
7015 new_node_pt->x(1) = x_tet[1];
7016 new_node_pt->x(2) = x_tet[2];
7017 }
7018 // Node already exists
7019 else
7020 {
7022 }
7023 }
7024
7025 // Mid brick-face node associated with face
7026 // spanned by mid-vertex nodes associated with tet edges 0 and 1
7027 {
7028 unsigned j = 4;
7029
7030 // Need new node?
7034
7036 if (old_node_pt == 0)
7037 {
7038 Node* new_node_pt = 0;
7039 if (face.is_boundary_face())
7040 {
7042 }
7043 else
7044 {
7046 }
7048 Node_pt.push_back(new_node_pt);
7049 Vector<double> s(3);
7053 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
7055 new_node_pt->x(0) = x_tet[0];
7056 new_node_pt->x(1) = x_tet[1];
7057 new_node_pt->x(2) = x_tet[2];
7058 }
7059 // Node already exists
7060 else
7061 {
7063 }
7064 }
7065
7066 // Top mid brick-face node -- only built by this element
7067 {
7068 unsigned j = 22;
7070 Node_pt.push_back(new_node_pt);
7071 Vector<double> s(3);
7075 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
7078 new_node_pt->x(0) = x_tet[0];
7079 new_node_pt->x(1) = x_tet[1];
7080 new_node_pt->x(2) = x_tet[2];
7081 }
7082
7083 // Right mid brick-face node -- only built by this element
7084 {
7085 unsigned j = 14;
7087 Node_pt.push_back(new_node_pt);
7088 Vector<double> s(3);
7092 dummy_q_el_pt[1]->interpolated_s_tet(s, s_tet);
7095 new_node_pt->x(0) = x_tet[0];
7096 new_node_pt->x(1) = x_tet[1];
7097 new_node_pt->x(2) = x_tet[2];
7098 }
7099
7100 // Back mid brick-face node copy from previous element
7101 {
7102 unsigned j = 16;
7103
7104 // Always copied
7106 }
7107 }
7108
7109 // Third brick element is centred at node 3 of tet:
7110 //-------------------------------------------------
7111 {
7112 // Assign coordinates of dummy element
7113 for (unsigned j = 0; j < 8; j++)
7114 {
7118 switch (j)
7119 {
7120 case 0:
7122 nod_pt->set_value(0, s_tet[0]);
7123 nod_pt->set_value(1, s_tet[1]);
7124 nod_pt->set_value(2, s_tet[2]);
7126 nod_pt->x(0) = x_tet[0];
7127 nod_pt->x(1) = x_tet[1];
7128 nod_pt->x(2) = x_tet[2];
7129 break;
7130 case 1:
7132 nod_pt->set_value(0, s_tet[0]);
7133 nod_pt->set_value(1, s_tet[1]);
7134 nod_pt->set_value(2, s_tet[2]);
7136 nod_pt->x(0) = x_tet[0];
7137 nod_pt->x(1) = x_tet[1];
7138 nod_pt->x(2) = x_tet[2];
7139 break;
7140 case 2:
7142 nod_pt->set_value(0, s_tet[0]);
7143 nod_pt->set_value(1, s_tet[1]);
7144 nod_pt->set_value(2, s_tet[2]);
7146 nod_pt->x(0) = x_tet[0];
7147 nod_pt->x(1) = x_tet[1];
7148 nod_pt->x(2) = x_tet[2];
7149 break;
7150 case 3:
7151 // label 13 in initial sketch: Mid face node on face
7152 // spanned by tet nodes 0,1,3
7153 s_tet[0] = 1.0 / 3.0;
7154 s_tet[1] = 1.0 / 3.0;
7155 s_tet[2] = 0.0;
7156 nod_pt->set_value(0, s_tet[0]);
7157 nod_pt->set_value(1, s_tet[1]);
7158 nod_pt->set_value(2, s_tet[2]);
7160 nod_pt->x(0) = x_tet[0];
7161 nod_pt->x(1) = x_tet[1];
7162 nod_pt->x(2) = x_tet[2];
7163 break;
7164 case 4:
7166 nod_pt->set_value(0, s_tet[0]);
7167 nod_pt->set_value(1, s_tet[1]);
7168 nod_pt->set_value(2, s_tet[2]);
7170 nod_pt->x(0) = x_tet[0];
7171 nod_pt->x(1) = x_tet[1];
7172 nod_pt->x(2) = x_tet[2];
7173 break;
7174 case 5:
7175 // label 12 in initial sketch: Mid face node on face
7176 // spanned by tet nodes 0,2,3
7177 s_tet[0] = 1.0 / 3.0;
7178 s_tet[1] = 0.0;
7179 s_tet[2] = 1.0 / 3.0;
7180 nod_pt->set_value(0, s_tet[0]);
7181 nod_pt->set_value(1, s_tet[1]);
7182 nod_pt->set_value(2, s_tet[2]);
7184 nod_pt->x(0) = x_tet[0];
7185 nod_pt->x(1) = x_tet[1];
7186 nod_pt->x(2) = x_tet[2];
7187 break;
7188 case 6:
7189 // label 10 in initial sketch: Mid face node on face
7190 // spanned by tet nodes 1,2,3
7191 s_tet[0] = 0.0;
7192 s_tet[1] = 1.0 / 3.0;
7193 s_tet[2] = 1.0 / 3.0;
7194 nod_pt->set_value(0, s_tet[0]);
7195 nod_pt->set_value(1, s_tet[1]);
7196 nod_pt->set_value(2, s_tet[2]);
7198 nod_pt->x(0) = x_tet[0];
7199 nod_pt->x(1) = x_tet[1];
7200 nod_pt->x(2) = x_tet[2];
7201 break;
7202 case 7:
7203 // label 14 in initial sketch: Centroid
7204 s_tet[0] = 0.25;
7205 s_tet[1] = 0.25;
7206 s_tet[2] = 0.25;
7207 nod_pt->set_value(0, s_tet[0]);
7208 nod_pt->set_value(1, s_tet[1]);
7209 nod_pt->set_value(2, s_tet[2]);
7211 nod_pt->x(0) = x_tet[0];
7212 nod_pt->x(1) = x_tet[1];
7213 nod_pt->x(2) = x_tet[2];
7214 break;
7215 }
7216 }
7217
7218 // Create actual second brick element
7219 FiniteElement* el_pt = new ELEMENT;
7221 Element_pt.push_back(el_pt);
7222
7223 TFace face0(
7225
7226 TFace face1(
7228
7229 TFace face2(
7231
7232 // Tet vertex nodes along edges emanating from node 0 in brick
7234 tet_edge_node[0].resize(2);
7235 tet_edge_node[0][0] = 6;
7236 tet_edge_node[0][1] = 0;
7237 tet_edge_node[1].resize(2);
7238 tet_edge_node[1][0] = 9;
7239 tet_edge_node[1][1] = 1;
7240 tet_edge_node[2].resize(2);
7241 tet_edge_node[2][0] = 8;
7242 tet_edge_node[2][1] = 2;
7243
7244 // Node number of tet vertex that node 0 in brick is centred on
7245 unsigned central_tet_vertex = 3;
7246
7247 Node* tet_node_pt = 0;
7248 Node* old_node_pt = 0;
7249
7250 // Corner node
7251 {
7252 unsigned j = 0;
7253
7254 // Need new node?
7257 if (old_node_pt == 0)
7258 {
7259 Node* new_node_pt = 0;
7260 if (tet_node_pt->is_on_boundary())
7261 {
7263 }
7264 else
7265 {
7267 }
7269 Node_pt.push_back(new_node_pt);
7270 Vector<double> s(3);
7274 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7276 new_node_pt->x(0) = x_tet[0];
7277 new_node_pt->x(1) = x_tet[1];
7278 new_node_pt->x(2) = x_tet[2];
7279 }
7280 // Node already exists
7281 else
7282 {
7284 }
7285 }
7286
7287 // Brick vertex node coindides with mid-edge node on tet edge 0
7288 {
7289 unsigned j = 2;
7290
7291 // Need new node?
7294 if (old_node_pt == 0)
7295 {
7296 Node* new_node_pt = 0;
7297 if (tet_node_pt->is_on_boundary())
7298 {
7300 }
7301 else
7302 {
7304 }
7306 Node_pt.push_back(new_node_pt);
7307 Vector<double> s(3);
7311 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7313 new_node_pt->x(0) = x_tet[0];
7314 new_node_pt->x(1) = x_tet[1];
7315 new_node_pt->x(2) = x_tet[2];
7316 }
7317 // Node already exists
7318 else
7319 {
7321 }
7322 }
7323
7324 // Brick vertex node coindides with mid vertex node of tet edge 1
7325 {
7326 unsigned j = 6;
7327
7328 // Need new node?
7331 if (old_node_pt == 0)
7332 {
7333 Node* new_node_pt = 0;
7334 if (tet_node_pt->is_on_boundary())
7335 {
7337 }
7338 else
7339 {
7341 }
7343 Node_pt.push_back(new_node_pt);
7344 Vector<double> s(3);
7348 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7350 new_node_pt->x(0) = x_tet[0];
7351 new_node_pt->x(1) = x_tet[1];
7352 new_node_pt->x(2) = x_tet[2];
7353 }
7354 // Node already exists
7355 else
7356 {
7358 }
7359 }
7360
7361 // Brick vertex node coindides with mid-vertex node of tet edge 2
7362 {
7363 unsigned j = 18;
7364
7365 // Need new node?
7368 if (old_node_pt == 0)
7369 {
7370 Node* new_node_pt = 0;
7371 if (tet_node_pt->is_on_boundary())
7372 {
7374 }
7375 else
7376 {
7378 }
7380 Node_pt.push_back(new_node_pt);
7381 Vector<double> s(3);
7385 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7387 new_node_pt->x(0) = x_tet[0];
7388 new_node_pt->x(1) = x_tet[1];
7389 new_node_pt->x(2) = x_tet[2];
7390 }
7391 // Node already exists
7392 else
7393 {
7395 }
7396 }
7397
7398 // Brick vertex node in the middle of tet face0
7399 {
7400 unsigned j = 20;
7401
7402 // Need new node?
7404 if (old_node_pt == 0)
7405 {
7406 Node* new_node_pt = 0;
7407 if (face0.is_boundary_face())
7408 {
7410 }
7411 else
7412 {
7414 }
7416 Node_pt.push_back(new_node_pt);
7417 Vector<double> s(3);
7421 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7423 new_node_pt->x(0) = x_tet[0];
7424 new_node_pt->x(1) = x_tet[1];
7425 new_node_pt->x(2) = x_tet[2];
7426 }
7427 // Node already exists
7428 else
7429 {
7431 }
7432 }
7433
7434 // Brick vertex node in the middle of tet face1
7435 {
7436 unsigned j = 24;
7437
7438 // Need new node?
7440 if (old_node_pt == 0)
7441 {
7442 Node* new_node_pt = 0;
7443 if (face1.is_boundary_face())
7444 {
7446 }
7447 else
7448 {
7450 }
7452 Node_pt.push_back(new_node_pt);
7453 Vector<double> s(3);
7457 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7459 new_node_pt->x(0) = x_tet[0];
7460 new_node_pt->x(1) = x_tet[1];
7461 new_node_pt->x(2) = x_tet[2];
7462 }
7463 // Node already exists
7464 else
7465 {
7467 }
7468 }
7469
7470 // Brick vertex node in the middle of tet face2
7471 {
7472 unsigned j = 8;
7473
7474 // Need new node?
7476 if (old_node_pt == 0)
7477 {
7478 Node* new_node_pt = 0;
7479 if (face2.is_boundary_face())
7480 {
7482 }
7483 else
7484 {
7486 }
7488 Node_pt.push_back(new_node_pt);
7489 Vector<double> s(3);
7493 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7495 new_node_pt->x(0) = x_tet[0];
7496 new_node_pt->x(1) = x_tet[1];
7497 new_node_pt->x(2) = x_tet[2];
7498 }
7499 // Node already exists
7500 else
7501 {
7503 }
7504 }
7505
7506 // Brick vertex node in centroid of tet. Only built for first element.
7507 // Enumerated "13" in initial sketch.
7508 {
7509 unsigned j = 26;
7510
7511 // Always copied
7513 }
7514
7515 // Internal brick node -- always built
7516 {
7517 unsigned j = 13;
7518
7519 // Always new
7520 {
7522 Node_pt.push_back(new_node_pt);
7523 Vector<double> s(3);
7527 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7529 new_node_pt->x(0) = x_tet[0];
7530 new_node_pt->x(1) = x_tet[1];
7531 new_node_pt->x(2) = x_tet[2];
7532 }
7533 }
7534
7535 // Brick edge node between brick nodes 0 and 2
7536 {
7537 unsigned j = 1;
7538
7539 // Need new node?
7540 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
7542 if (old_node_pt == 0)
7543 {
7544 Node* new_node_pt = 0;
7545 if (edge.is_boundary_edge())
7546 {
7548 }
7549 else
7550 {
7552 }
7554 Node_pt.push_back(new_node_pt);
7555 Vector<double> s(3);
7559 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7561 new_node_pt->x(0) = x_tet[0];
7562 new_node_pt->x(1) = x_tet[1];
7563 new_node_pt->x(2) = x_tet[2];
7564 }
7565 // Node already exists
7566 else
7567 {
7569 }
7570 }
7571
7572 // Brick edge node between brick nodes 0 and 6
7573 {
7574 unsigned j = 3;
7575
7576 // Need new node?
7577 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
7579 if (old_node_pt == 0)
7580 {
7581 Node* new_node_pt = 0;
7582 if (edge.is_boundary_edge())
7583 {
7585 }
7586 else
7587 {
7589 }
7591 Node_pt.push_back(new_node_pt);
7592 Vector<double> s(3);
7596 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7598 new_node_pt->x(0) = x_tet[0];
7599 new_node_pt->x(1) = x_tet[1];
7600 new_node_pt->x(2) = x_tet[2];
7601 }
7602 // Node already exists
7603 else
7604 {
7606 }
7607 }
7608
7609 // Brick edge node between brick nodes 2 and 8
7610 {
7611 unsigned j = 5;
7612
7613 // Need new node?
7614 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
7616 if (old_node_pt == 0)
7617 {
7618 Node* new_node_pt = 0;
7619 if (edge.is_boundary_edge())
7620 {
7622 }
7623 else
7624 {
7626 }
7628 Node_pt.push_back(new_node_pt);
7629 Vector<double> s(3);
7633 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7635 new_node_pt->x(0) = x_tet[0];
7636 new_node_pt->x(1) = x_tet[1];
7637 new_node_pt->x(2) = x_tet[2];
7638 }
7639 // Node already exists
7640 else
7641 {
7643 }
7644 }
7645
7646 // Brick edge node between brick nodes 6 and 8
7647 {
7648 unsigned j = 7;
7649
7650 // Need new node?
7651 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
7653 if (old_node_pt == 0)
7654 {
7655 Node* new_node_pt = 0;
7656 if (edge.is_boundary_edge())
7657 {
7659 }
7660 else
7661 {
7663 }
7665 Node_pt.push_back(new_node_pt);
7666 Vector<double> s(3);
7670 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7672 new_node_pt->x(0) = x_tet[0];
7673 new_node_pt->x(1) = x_tet[1];
7674 new_node_pt->x(2) = x_tet[2];
7675 }
7676 // Node already exists
7677 else
7678 {
7680 }
7681 }
7682
7683 // Brick edge node between brick nodes 18 and 20
7684 {
7685 unsigned j = 19;
7686
7687 // Need new node?
7688 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
7690 if (old_node_pt == 0)
7691 {
7692 Node* new_node_pt = 0;
7693 if (edge.is_boundary_edge())
7694 {
7696 }
7697 else
7698 {
7700 }
7702 Node_pt.push_back(new_node_pt);
7703 Vector<double> s(3);
7707 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7709 new_node_pt->x(0) = x_tet[0];
7710 new_node_pt->x(1) = x_tet[1];
7711 new_node_pt->x(2) = x_tet[2];
7712 }
7713 // Node already exists
7714 else
7715 {
7717 }
7718 }
7719
7720 // Brick edge node between brick nodes 18 and 24
7721 {
7722 unsigned j = 21;
7723
7724 // Need new node?
7725 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
7727 if (old_node_pt == 0)
7728 {
7729 Node* new_node_pt = 0;
7730 if (edge.is_boundary_edge())
7731 {
7733 }
7734 else
7735 {
7737 }
7739 Node_pt.push_back(new_node_pt);
7740 Vector<double> s(3);
7744 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7746 new_node_pt->x(0) = x_tet[0];
7747 new_node_pt->x(1) = x_tet[1];
7748 new_node_pt->x(2) = x_tet[2];
7749 }
7750 // Node already exists
7751 else
7752 {
7754 }
7755 }
7756
7757 // Brick edge node between brick nodes 20 and 26
7758 {
7759 unsigned j = 23;
7760
7761 // Need new node?
7762 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
7764 if (old_node_pt == 0)
7765 {
7766 Node* new_node_pt = 0;
7767 if (edge.is_boundary_edge())
7768 {
7770 }
7771 else
7772 {
7774 }
7776 Node_pt.push_back(new_node_pt);
7777 Vector<double> s(3);
7781 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7783 new_node_pt->x(0) = x_tet[0];
7784 new_node_pt->x(1) = x_tet[1];
7785 new_node_pt->x(2) = x_tet[2];
7786 }
7787 // Node already exists
7788 else
7789 {
7791 }
7792 }
7793
7794 // Brick edge node between brick nodes 24 and 26
7795 {
7796 unsigned j = 25;
7797
7798 // Need new node?
7799 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
7801 if (old_node_pt == 0)
7802 {
7803 Node* new_node_pt = 0;
7804 if (edge.is_boundary_edge())
7805 {
7807 }
7808 else
7809 {
7811 }
7813 Node_pt.push_back(new_node_pt);
7814 Vector<double> s(3);
7818 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7820 new_node_pt->x(0) = x_tet[0];
7821 new_node_pt->x(1) = x_tet[1];
7822 new_node_pt->x(2) = x_tet[2];
7823 }
7824 // Node already exists
7825 else
7826 {
7828 }
7829 }
7830
7831 // Brick edge node between brick nodes 0 and 18
7832 {
7833 unsigned j = 9;
7834
7835 // Need new node?
7836 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
7838 if (old_node_pt == 0)
7839 {
7840 Node* new_node_pt = 0;
7841 if (edge.is_boundary_edge())
7842 {
7844 }
7845 else
7846 {
7848 }
7850 Node_pt.push_back(new_node_pt);
7851 Vector<double> s(3);
7855 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7857 new_node_pt->x(0) = x_tet[0];
7858 new_node_pt->x(1) = x_tet[1];
7859 new_node_pt->x(2) = x_tet[2];
7860 }
7861 // Node already exists
7862 else
7863 {
7865 }
7866 }
7867
7868 // Brick edge node between brick nodes 2 and 20
7869 {
7870 unsigned j = 11;
7871
7872 // Need new node?
7873 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
7875 if (old_node_pt == 0)
7876 {
7877 Node* new_node_pt = 0;
7878 if (edge.is_boundary_edge())
7879 {
7881 }
7882 else
7883 {
7885 }
7887 Node_pt.push_back(new_node_pt);
7888 Vector<double> s(3);
7892 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7894 new_node_pt->x(0) = x_tet[0];
7895 new_node_pt->x(1) = x_tet[1];
7896 new_node_pt->x(2) = x_tet[2];
7897 }
7898 // Node already exists
7899 else
7900 {
7902 }
7903 }
7904
7905 // Brick edge node between brick nodes 6 and 24
7906 {
7907 unsigned j = 15;
7908
7909 // Need new node?
7910 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
7912 if (old_node_pt == 0)
7913 {
7914 Node* new_node_pt = 0;
7915 if (edge.is_boundary_edge())
7916 {
7918 }
7919 else
7920 {
7922 }
7924 Node_pt.push_back(new_node_pt);
7925 Vector<double> s(3);
7929 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7931 new_node_pt->x(0) = x_tet[0];
7932 new_node_pt->x(1) = x_tet[1];
7933 new_node_pt->x(2) = x_tet[2];
7934 }
7935 // Node already exists
7936 else
7937 {
7939 }
7940 }
7941
7942 // Brick edge node between brick nodes 8 and 26
7943 {
7944 unsigned j = 17;
7945
7946 // Need new node?
7947 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
7949 if (old_node_pt == 0)
7950 {
7951 Node* new_node_pt = 0;
7952 if (edge.is_boundary_edge())
7953 {
7955 }
7956 else
7957 {
7959 }
7961 Node_pt.push_back(new_node_pt);
7962 Vector<double> s(3);
7966 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
7968 new_node_pt->x(0) = x_tet[0];
7969 new_node_pt->x(1) = x_tet[1];
7970 new_node_pt->x(2) = x_tet[2];
7971 }
7972 // Node already exists
7973 else
7974 {
7976 }
7977 }
7978
7979 // Mid brick-face node associated with face
7980 // spanned by mid-vertex nodes associated with tet edges 0 and 2
7981 {
7982 unsigned j = 10;
7983
7984 // Need new node?
7988
7990 if (old_node_pt == 0)
7991 {
7992 Node* new_node_pt = 0;
7993 if (face.is_boundary_face())
7994 {
7996 }
7997 else
7998 {
8000 }
8002 Node_pt.push_back(new_node_pt);
8003 Vector<double> s(3);
8007 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
8009 new_node_pt->x(0) = x_tet[0];
8010 new_node_pt->x(1) = x_tet[1];
8011 new_node_pt->x(2) = x_tet[2];
8012 }
8013 // Node already exists
8014 else
8015 {
8017 }
8018 }
8019
8020 // Mid brick-face node associated with face
8021 // spanned by mid-vertex nodes associated with tet edges 1 and 2
8022 {
8023 unsigned j = 12;
8024
8025 // Need new node?
8030 if (old_node_pt == 0)
8031 {
8032 Node* new_node_pt = 0;
8033 if (face.is_boundary_face())
8034 {
8036 }
8037 else
8038 {
8040 }
8042 Node_pt.push_back(new_node_pt);
8043 Vector<double> s(3);
8047 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
8049 new_node_pt->x(0) = x_tet[0];
8050 new_node_pt->x(1) = x_tet[1];
8051 new_node_pt->x(2) = x_tet[2];
8052 }
8053 // Node already exists
8054 else
8055 {
8057 }
8058 }
8059
8060 // Mid brick-face node associated with face
8061 // spanned by mid-vertex nodes associated with tet edges 0 and 1
8062 {
8063 unsigned j = 4;
8064
8065 // Need new node?
8069
8071 if (old_node_pt == 0)
8072 {
8073 Node* new_node_pt = 0;
8074 if (face.is_boundary_face())
8075 {
8077 }
8078 else
8079 {
8081 }
8083 Node_pt.push_back(new_node_pt);
8084 Vector<double> s(3);
8088 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
8090 new_node_pt->x(0) = x_tet[0];
8091 new_node_pt->x(1) = x_tet[1];
8092 new_node_pt->x(2) = x_tet[2];
8093 }
8094 // Node already exists
8095 else
8096 {
8098 }
8099 }
8100
8101 // Top mid brick-face node -- only built by this element
8102 {
8103 unsigned j = 22;
8105 Node_pt.push_back(new_node_pt);
8106 Vector<double> s(3);
8110 dummy_q_el_pt[2]->interpolated_s_tet(s, s_tet);
8113 new_node_pt->x(0) = x_tet[0];
8114 new_node_pt->x(1) = x_tet[1];
8115 new_node_pt->x(2) = x_tet[2];
8116 }
8117
8118 // Right mid brick-face node copy from first element
8119 {
8120 unsigned j = 14;
8121
8122 // Always copied
8124 }
8125
8126 // Back mid brick-face node copy from previous element
8127 {
8128 unsigned j = 16;
8129
8130 // Always copied
8132 }
8133 }
8134
8135 // Fourth brick element is centred at node 2 of tet:
8136 //--------------------------------------------------
8137 {
8138 // Assign coordinates of dummy element
8139 for (unsigned j = 0; j < 8; j++)
8140 {
8144 switch (j)
8145 {
8146 case 0:
8148 nod_pt->set_value(0, s_tet[0]);
8149 nod_pt->set_value(1, s_tet[1]);
8150 nod_pt->set_value(2, s_tet[2]);
8152 nod_pt->x(0) = x_tet[0];
8153 nod_pt->x(1) = x_tet[1];
8154 nod_pt->x(2) = x_tet[2];
8155 break;
8156 case 1:
8158 nod_pt->set_value(0, s_tet[0]);
8159 nod_pt->set_value(1, s_tet[1]);
8160 nod_pt->set_value(2, s_tet[2]);
8162 nod_pt->x(0) = x_tet[0];
8163 nod_pt->x(1) = x_tet[1];
8164 nod_pt->x(2) = x_tet[2];
8165 break;
8166 case 2:
8168 nod_pt->set_value(0, s_tet[0]);
8169 nod_pt->set_value(1, s_tet[1]);
8170 nod_pt->set_value(2, s_tet[2]);
8172 nod_pt->x(0) = x_tet[0];
8173 nod_pt->x(1) = x_tet[1];
8174 nod_pt->x(2) = x_tet[2];
8175 break;
8176 case 3:
8177 // label 11 in initial sketch: Mid face node on face
8178 // spanned by tet nodes 0,1,2
8179 s_tet[0] = 1.0 / 3.0;
8180 s_tet[1] = 1.0 / 3.0;
8181 s_tet[2] = 1.0 / 3.0;
8182 nod_pt->set_value(0, s_tet[0]);
8183 nod_pt->set_value(1, s_tet[1]);
8184 nod_pt->set_value(2, s_tet[2]);
8186 nod_pt->x(0) = x_tet[0];
8187 nod_pt->x(1) = x_tet[1];
8188 nod_pt->x(2) = x_tet[2];
8189 break;
8190 case 4:
8192 nod_pt->set_value(0, s_tet[0]);
8193 nod_pt->set_value(1, s_tet[1]);
8194 nod_pt->set_value(2, s_tet[2]);
8196 nod_pt->x(0) = x_tet[0];
8197 nod_pt->x(1) = x_tet[1];
8198 nod_pt->x(2) = x_tet[2];
8199 break;
8200 case 5:
8201 // label 10 in initial sketch: Mid face node on face
8202 // spanned by tet nodes 0,2,3
8203 s_tet[0] = 0.0;
8204 s_tet[1] = 1.0 / 3.0;
8205 s_tet[2] = 1.0 / 3.0;
8206 nod_pt->set_value(0, s_tet[0]);
8207 nod_pt->set_value(1, s_tet[1]);
8208 nod_pt->set_value(2, s_tet[2]);
8210 nod_pt->x(0) = x_tet[0];
8211 nod_pt->x(1) = x_tet[1];
8212 nod_pt->x(2) = x_tet[2];
8213 break;
8214 case 6:
8215 // label 12 in initial sketch: Mid face node on face
8216 // spanned by tet nodes 0,2,3
8217 s_tet[0] = 1.0 / 3.0;
8218 s_tet[1] = 0.0;
8219 s_tet[2] = 1.0 / 3.0;
8220 nod_pt->set_value(0, s_tet[0]);
8221 nod_pt->set_value(1, s_tet[1]);
8222 nod_pt->set_value(2, s_tet[2]);
8224 nod_pt->x(0) = x_tet[0];
8225 nod_pt->x(1) = x_tet[1];
8226 nod_pt->x(2) = x_tet[2];
8227 break;
8228 case 7:
8229 // label 14 in initial sketch: Centroid
8230 s_tet[0] = 0.25;
8231 s_tet[1] = 0.25;
8232 s_tet[2] = 0.25;
8233 nod_pt->set_value(0, s_tet[0]);
8234 nod_pt->set_value(1, s_tet[1]);
8235 nod_pt->set_value(2, s_tet[2]);
8237 nod_pt->x(0) = x_tet[0];
8238 nod_pt->x(1) = x_tet[1];
8239 nod_pt->x(2) = x_tet[2];
8240 break;
8241 }
8242 }
8243
8244 // Create actual third brick element
8245 FiniteElement* el_pt = new ELEMENT;
8247 Element_pt.push_back(el_pt);
8248
8249 TFace face0(
8251
8252 TFace face1(
8254
8255 TFace face2(
8257
8258 // Tet vertex nodes along edges emanating from node 0 in brick
8260 tet_edge_node[0].resize(2);
8261 tet_edge_node[0][0] = 7;
8262 tet_edge_node[0][1] = 1;
8263 tet_edge_node[1].resize(2);
8264 tet_edge_node[1][0] = 5;
8265 tet_edge_node[1][1] = 0;
8266 tet_edge_node[2].resize(2);
8267 tet_edge_node[2][0] = 8;
8268 tet_edge_node[2][1] = 3;
8269
8270 // Node number of tet vertex that node 0 in brick is centred on
8271 unsigned central_tet_vertex = 2;
8272
8273 Node* tet_node_pt = 0;
8274 Node* old_node_pt = 0;
8275
8276 // Corner node
8277 {
8278 unsigned j = 0;
8279
8280 // Need new node?
8283 if (old_node_pt == 0)
8284 {
8285 Node* new_node_pt = 0;
8286 if (tet_node_pt->is_on_boundary())
8287 {
8289 }
8290 else
8291 {
8293 }
8295 Node_pt.push_back(new_node_pt);
8296 Vector<double> s(3);
8300 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8302 new_node_pt->x(0) = x_tet[0];
8303 new_node_pt->x(1) = x_tet[1];
8304 new_node_pt->x(2) = x_tet[2];
8305 }
8306 // Node already exists
8307 else
8308 {
8310 }
8311 }
8312
8313 // Brick vertex node coindides with mid-edge node on tet edge 0
8314 {
8315 unsigned j = 2;
8316
8317 // Need new node?
8320 if (old_node_pt == 0)
8321 {
8322 Node* new_node_pt = 0;
8323 if (tet_node_pt->is_on_boundary())
8324 {
8326 }
8327 else
8328 {
8330 }
8332 Node_pt.push_back(new_node_pt);
8333 Vector<double> s(3);
8337 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8339 new_node_pt->x(0) = x_tet[0];
8340 new_node_pt->x(1) = x_tet[1];
8341 new_node_pt->x(2) = x_tet[2];
8342 }
8343 // Node already exists
8344 else
8345 {
8347 }
8348 }
8349
8350 // Brick vertex node coindides with mid vertex node of tet edge 1
8351 {
8352 unsigned j = 6;
8353
8354 // Need new node?
8357 if (old_node_pt == 0)
8358 {
8359 Node* new_node_pt = 0;
8360 if (tet_node_pt->is_on_boundary())
8361 {
8363 }
8364 else
8365 {
8367 }
8369 Node_pt.push_back(new_node_pt);
8370 Vector<double> s(3);
8374 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8376 new_node_pt->x(0) = x_tet[0];
8377 new_node_pt->x(1) = x_tet[1];
8378 new_node_pt->x(2) = x_tet[2];
8379 }
8380 // Node already exists
8381 else
8382 {
8384 }
8385 }
8386
8387 // Brick vertex node coindides with mid-vertex node of tet edge 2
8388 {
8389 unsigned j = 18;
8390
8391 // Need new node?
8394 if (old_node_pt == 0)
8395 {
8396 Node* new_node_pt = 0;
8397 if (tet_node_pt->is_on_boundary())
8398 {
8400 }
8401 else
8402 {
8404 }
8406 Node_pt.push_back(new_node_pt);
8407 Vector<double> s(3);
8411 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8413 new_node_pt->x(0) = x_tet[0];
8414 new_node_pt->x(1) = x_tet[1];
8415 new_node_pt->x(2) = x_tet[2];
8416 }
8417 // Node already exists
8418 else
8419 {
8421 }
8422 }
8423
8424 // Brick vertex node in the middle of tet face0
8425 {
8426 unsigned j = 20;
8427
8428 // Need new node?
8430 if (old_node_pt == 0)
8431 {
8432 Node* new_node_pt = 0;
8433 if (face0.is_boundary_face())
8434 {
8436 }
8437 else
8438 {
8440 }
8442 Node_pt.push_back(new_node_pt);
8443 Vector<double> s(3);
8447 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8449 new_node_pt->x(0) = x_tet[0];
8450 new_node_pt->x(1) = x_tet[1];
8451 new_node_pt->x(2) = x_tet[2];
8452 }
8453 // Node already exists
8454 else
8455 {
8457 }
8458 }
8459
8460 // Brick vertex node in the middle of tet face1
8461 {
8462 unsigned j = 24;
8463
8464 // Need new node?
8466 if (old_node_pt == 0)
8467 {
8468 Node* new_node_pt = 0;
8469 if (face1.is_boundary_face())
8470 {
8472 }
8473 else
8474 {
8476 }
8478 Node_pt.push_back(new_node_pt);
8479 Vector<double> s(3);
8483 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8485 new_node_pt->x(0) = x_tet[0];
8486 new_node_pt->x(1) = x_tet[1];
8487 new_node_pt->x(2) = x_tet[2];
8488 }
8489 // Node already exists
8490 else
8491 {
8493 }
8494 }
8495
8496 // Brick vertex node in the middle of tet face2
8497 {
8498 unsigned j = 8;
8499
8500 // Need new node?
8502 if (old_node_pt == 0)
8503 {
8504 Node* new_node_pt = 0;
8505 if (face2.is_boundary_face())
8506 {
8508 }
8509 else
8510 {
8512 }
8514 Node_pt.push_back(new_node_pt);
8515 Vector<double> s(3);
8519 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8521 new_node_pt->x(0) = x_tet[0];
8522 new_node_pt->x(1) = x_tet[1];
8523 new_node_pt->x(2) = x_tet[2];
8524 }
8525 // Node already exists
8526 else
8527 {
8529 }
8530 }
8531
8532 // Brick vertex node in centroid of tet. Only built for first element.
8533 // Enumerated "13" in initial sketch.
8534 {
8535 unsigned j = 26;
8536
8537 // Always copied
8539 }
8540
8541 // Internal brick node -- always built
8542 {
8543 unsigned j = 13;
8544
8545 // Always new
8546 {
8548 Node_pt.push_back(new_node_pt);
8549 Vector<double> s(3);
8553 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8555 new_node_pt->x(0) = x_tet[0];
8556 new_node_pt->x(1) = x_tet[1];
8557 new_node_pt->x(2) = x_tet[2];
8558 }
8559 }
8560
8561 // Brick edge node between brick nodes 0 and 2
8562 {
8563 unsigned j = 1;
8564
8565 // Need new node?
8566 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
8568 if (old_node_pt == 0)
8569 {
8570 Node* new_node_pt = 0;
8571 if (edge.is_boundary_edge())
8572 {
8574 }
8575 else
8576 {
8578 }
8580 Node_pt.push_back(new_node_pt);
8581 Vector<double> s(3);
8585 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8587 new_node_pt->x(0) = x_tet[0];
8588 new_node_pt->x(1) = x_tet[1];
8589 new_node_pt->x(2) = x_tet[2];
8590 }
8591 // Node already exists
8592 else
8593 {
8595 }
8596 }
8597
8598 // Brick edge node between brick nodes 0 and 6
8599 {
8600 unsigned j = 3;
8601
8602 // Need new node?
8603 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
8605 if (old_node_pt == 0)
8606 {
8607 Node* new_node_pt = 0;
8608 if (edge.is_boundary_edge())
8609 {
8611 }
8612 else
8613 {
8615 }
8617 Node_pt.push_back(new_node_pt);
8618 Vector<double> s(3);
8622 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8624 new_node_pt->x(0) = x_tet[0];
8625 new_node_pt->x(1) = x_tet[1];
8626 new_node_pt->x(2) = x_tet[2];
8627 }
8628 // Node already exists
8629 else
8630 {
8632 }
8633 }
8634
8635 // Brick edge node between brick nodes 2 and 8
8636 {
8637 unsigned j = 5;
8638
8639 // Need new node?
8640 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
8642 if (old_node_pt == 0)
8643 {
8644 Node* new_node_pt = 0;
8645 if (edge.is_boundary_edge())
8646 {
8648 }
8649 else
8650 {
8652 }
8654 Node_pt.push_back(new_node_pt);
8655 Vector<double> s(3);
8659 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8661 new_node_pt->x(0) = x_tet[0];
8662 new_node_pt->x(1) = x_tet[1];
8663 new_node_pt->x(2) = x_tet[2];
8664 }
8665 // Node already exists
8666 else
8667 {
8669 }
8670 }
8671
8672 // Brick edge node between brick nodes 6 and 8
8673 {
8674 unsigned j = 7;
8675
8676 // Need new node?
8677 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
8679 if (old_node_pt == 0)
8680 {
8681 Node* new_node_pt = 0;
8682 if (edge.is_boundary_edge())
8683 {
8685 }
8686 else
8687 {
8689 }
8691 Node_pt.push_back(new_node_pt);
8692 Vector<double> s(3);
8696 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8698 new_node_pt->x(0) = x_tet[0];
8699 new_node_pt->x(1) = x_tet[1];
8700 new_node_pt->x(2) = x_tet[2];
8701 }
8702 // Node already exists
8703 else
8704 {
8706 }
8707 }
8708
8709 // Brick edge node between brick nodes 18 and 20
8710 {
8711 unsigned j = 19;
8712
8713 // Need new node?
8714 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
8716 if (old_node_pt == 0)
8717 {
8718 Node* new_node_pt = 0;
8719 if (edge.is_boundary_edge())
8720 {
8722 }
8723 else
8724 {
8726 }
8728 Node_pt.push_back(new_node_pt);
8729 Vector<double> s(3);
8733 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8735 new_node_pt->x(0) = x_tet[0];
8736 new_node_pt->x(1) = x_tet[1];
8737 new_node_pt->x(2) = x_tet[2];
8738 }
8739 // Node already exists
8740 else
8741 {
8743 }
8744 }
8745
8746 // Brick edge node between brick nodes 18 and 24
8747 {
8748 unsigned j = 21;
8749
8750 // Need new node?
8751 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
8753 if (old_node_pt == 0)
8754 {
8755 Node* new_node_pt = 0;
8756 if (edge.is_boundary_edge())
8757 {
8759 }
8760 else
8761 {
8763 }
8765 Node_pt.push_back(new_node_pt);
8766 Vector<double> s(3);
8770 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8772 new_node_pt->x(0) = x_tet[0];
8773 new_node_pt->x(1) = x_tet[1];
8774 new_node_pt->x(2) = x_tet[2];
8775 }
8776 // Node already exists
8777 else
8778 {
8780 }
8781 }
8782
8783 // Brick edge node between brick nodes 20 and 26
8784 {
8785 unsigned j = 23;
8786
8787 // Need new node?
8788 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
8790 if (old_node_pt == 0)
8791 {
8792 Node* new_node_pt = 0;
8793 if (edge.is_boundary_edge())
8794 {
8796 }
8797 else
8798 {
8800 }
8802 Node_pt.push_back(new_node_pt);
8803 Vector<double> s(3);
8807 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8809 new_node_pt->x(0) = x_tet[0];
8810 new_node_pt->x(1) = x_tet[1];
8811 new_node_pt->x(2) = x_tet[2];
8812 }
8813 // Node already exists
8814 else
8815 {
8817 }
8818 }
8819
8820 // Brick edge node between brick nodes 24 and 26
8821 {
8822 unsigned j = 25;
8823
8824 // Need new node?
8825 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
8827 if (old_node_pt == 0)
8828 {
8829 Node* new_node_pt = 0;
8830 if (edge.is_boundary_edge())
8831 {
8833 }
8834 else
8835 {
8837 }
8839 Node_pt.push_back(new_node_pt);
8840 Vector<double> s(3);
8844 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8846 new_node_pt->x(0) = x_tet[0];
8847 new_node_pt->x(1) = x_tet[1];
8848 new_node_pt->x(2) = x_tet[2];
8849 }
8850 // Node already exists
8851 else
8852 {
8854 }
8855 }
8856
8857 // Brick edge node between brick nodes 0 and 18
8858 {
8859 unsigned j = 9;
8860
8861 // Need new node?
8862 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
8864 if (old_node_pt == 0)
8865 {
8866 Node* new_node_pt = 0;
8867 if (edge.is_boundary_edge())
8868 {
8870 }
8871 else
8872 {
8874 }
8876 Node_pt.push_back(new_node_pt);
8877 Vector<double> s(3);
8881 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8883 new_node_pt->x(0) = x_tet[0];
8884 new_node_pt->x(1) = x_tet[1];
8885 new_node_pt->x(2) = x_tet[2];
8886 }
8887 // Node already exists
8888 else
8889 {
8891 }
8892 }
8893
8894 // Brick edge node between brick nodes 2 and 20
8895 {
8896 unsigned j = 11;
8897
8898 // Need new node?
8899 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
8901 if (old_node_pt == 0)
8902 {
8903 Node* new_node_pt = 0;
8904 if (edge.is_boundary_edge())
8905 {
8907 }
8908 else
8909 {
8911 }
8913 Node_pt.push_back(new_node_pt);
8914 Vector<double> s(3);
8918 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8920 new_node_pt->x(0) = x_tet[0];
8921 new_node_pt->x(1) = x_tet[1];
8922 new_node_pt->x(2) = x_tet[2];
8923 }
8924 // Node already exists
8925 else
8926 {
8928 }
8929 }
8930
8931 // Brick edge node between brick nodes 6 and 24
8932 {
8933 unsigned j = 15;
8934
8935 // Need new node?
8936 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
8938 if (old_node_pt == 0)
8939 {
8940 Node* new_node_pt = 0;
8941 if (edge.is_boundary_edge())
8942 {
8944 }
8945 else
8946 {
8948 }
8950 Node_pt.push_back(new_node_pt);
8951 Vector<double> s(3);
8955 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8957 new_node_pt->x(0) = x_tet[0];
8958 new_node_pt->x(1) = x_tet[1];
8959 new_node_pt->x(2) = x_tet[2];
8960 }
8961 // Node already exists
8962 else
8963 {
8965 }
8966 }
8967
8968 // Brick edge node between brick nodes 8 and 26
8969 {
8970 unsigned j = 17;
8971
8972 // Need new node?
8973 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
8975 if (old_node_pt == 0)
8976 {
8977 Node* new_node_pt = 0;
8978 if (edge.is_boundary_edge())
8979 {
8981 }
8982 else
8983 {
8985 }
8987 Node_pt.push_back(new_node_pt);
8988 Vector<double> s(3);
8992 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
8994 new_node_pt->x(0) = x_tet[0];
8995 new_node_pt->x(1) = x_tet[1];
8996 new_node_pt->x(2) = x_tet[2];
8997 }
8998 // Node already exists
8999 else
9000 {
9002 }
9003 }
9004
9005 // Mid brick-face node associated with face
9006 // spanned by mid-vertex nodes associated with tet edges 0 and 2
9007 {
9008 unsigned j = 10;
9009
9010 // Need new node?
9014
9016 if (old_node_pt == 0)
9017 {
9018 Node* new_node_pt = 0;
9019 if (face.is_boundary_face())
9020 {
9022 }
9023 else
9024 {
9026 }
9028 Node_pt.push_back(new_node_pt);
9029 Vector<double> s(3);
9033 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
9035 new_node_pt->x(0) = x_tet[0];
9036 new_node_pt->x(1) = x_tet[1];
9037 new_node_pt->x(2) = x_tet[2];
9038 }
9039 // Node already exists
9040 else
9041 {
9043 }
9044 }
9045
9046 // Mid brick-face node associated with face
9047 // spanned by mid-vertex nodes associated with tet edges 1 and 2
9048 {
9049 unsigned j = 12;
9050
9051 // Need new node?
9055
9057 if (old_node_pt == 0)
9058 {
9059 Node* new_node_pt = 0;
9060 if (face.is_boundary_face())
9061 {
9063 }
9064 else
9065 {
9067 }
9069 Node_pt.push_back(new_node_pt);
9070 Vector<double> s(3);
9074 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
9076 new_node_pt->x(0) = x_tet[0];
9077 new_node_pt->x(1) = x_tet[1];
9078 new_node_pt->x(2) = x_tet[2];
9079 }
9080 // Node already exists
9081 else
9082 {
9084 }
9085 }
9086
9087 // Mid brick-face node associated with face
9088 // spanned by mid-vertex nodes associated with tet edges 0 and 1
9089 {
9090 unsigned j = 4;
9091
9092 // Need new node?
9096
9098 if (old_node_pt == 0)
9099 {
9100 Node* new_node_pt = 0;
9101 if (face.is_boundary_face())
9102 {
9104 }
9105 else
9106 {
9108 }
9110 Node_pt.push_back(new_node_pt);
9111 Vector<double> s(3);
9115 dummy_q_el_pt[3]->interpolated_s_tet(s, s_tet);
9117 new_node_pt->x(0) = x_tet[0];
9118 new_node_pt->x(1) = x_tet[1];
9119 new_node_pt->x(2) = x_tet[2];
9120 }
9121 // Node already exists
9122 else
9123 {
9125 }
9126 }
9127
9128 // Top mid brick-face node copy from top of second element
9129 {
9130 unsigned j = 22;
9131
9132 // Always copied
9134 }
9135
9136 // Right mid brick-face node copy from top of first element
9137 {
9138 unsigned j = 14;
9139
9140 // Always copied
9142 }
9143
9144 // Back mid brick-face node copy from top of zeroth element
9145 {
9146 unsigned j = 16;
9147
9148 // Always copied
9150 }
9151 }
9152
9153 // Check if the four faces of the tet are on a boundary.
9154 // If they are, add the nodes in the brick mesh to those
9155 // boundaries.
9156
9157 // Loop over faces of tet (in its own face index enumeration)
9158 for (int face_index = 0; face_index < 4; face_index++)
9159 {
9160 // Identify face and coordinates in it
9161 TFace* face_pt = 0;
9162 switch (face_index)
9163 {
9164 case 0:
9165 // Face 0: s[0]=0; opposite node 0
9166 face_pt = new TFace(tet_el_pt->node_pt(1),
9167 tet_el_pt->node_pt(2),
9168 tet_el_pt->node_pt(3));
9169 break;
9170
9171 case 1:
9172 // Face 1: s[1]=0; opposite node 1
9173 face_pt = new TFace(tet_el_pt->node_pt(0),
9174 tet_el_pt->node_pt(2),
9175 tet_el_pt->node_pt(3));
9176
9177 break;
9178
9179 case 2:
9180 // Face 2: s[2]=0; opposite node 2
9181 face_pt = new TFace(tet_el_pt->node_pt(0),
9182 tet_el_pt->node_pt(1),
9183 tet_el_pt->node_pt(3));
9184
9185 break;
9186
9187 case 3:
9188 // Face 3: s[0]+s[1]+s[2]=1; opposite node 3
9189 face_pt = new TFace(tet_el_pt->node_pt(0),
9190 tet_el_pt->node_pt(1),
9191 tet_el_pt->node_pt(2));
9192 break;
9193 }
9194
9195 if (face_pt->is_boundary_face())
9196 {
9197 std::set<unsigned> bnd;
9198 std::set<unsigned>* bnd_pt = &bnd;
9199 face_pt->get_boundaries_pt(bnd_pt);
9200
9201#ifdef PARANOID
9202 if ((*bnd_pt).size() > 1)
9203 {
9204 std::ostringstream error_stream;
9205 error_stream << "TFace should only be on one boundary.\n";
9206 throw OomphLibError(error_stream.str(),
9209 }
9210#endif
9211
9212 if ((*bnd_pt).size() == 1)
9213 {
9214 // Create new face element
9217 {
9218#ifdef PARANOID
9219 if (dynamic_cast<SolidTElement<3, 3>*>(tet_el_pt) == 0)
9220 {
9221 std::ostringstream error_stream;
9223 << "Tet-element cannot be cast to SolidTElement<3,3>.\n"
9224 << "BrickFromTetMesh can only be built from\n"
9225 << "mesh containing quadratic tets.\n"
9226 << std::endl;
9227 throw OomphLibError(error_stream.str(),
9230 }
9231#endif
9232 // Build the face element
9234 tet_el_pt, face_index);
9235 }
9236 else
9237 {
9238#ifdef PARANOID
9239 if (dynamic_cast<TElement<3, 3>*>(tet_el_pt) == 0)
9240 {
9241 std::ostringstream error_stream;
9242 error_stream << "Tet-element cannot be cast to TElement<3,3>.\n"
9243 << "BrickFromTetMesh can only be built from\n"
9244 << "mesh containing quadratic tets.\n"
9245 << std::endl;
9246 throw OomphLibError(error_stream.str(),
9249 }
9250#endif
9251
9252 // Build the face element
9253 face_el_pt =
9255 }
9256
9257 // Specify boundary id in bulk mesh (needed to extract
9258 // boundary coordinate)
9259 unsigned b = (*(*bnd_pt).begin());
9260 set_boundary_coordinate_exists(b);
9261 face_el_pt->set_boundary_number_in_bulk_mesh(b);
9262
9263 // Now set up the brick nodes on this face, enumerated as
9264 // in s_face
9266
9267 switch (face_index)
9268 {
9269 case 0:
9273
9277
9281
9285
9287
9290
9293
9296 break;
9297
9298 case 1:
9302
9306
9310
9314
9316
9319
9322
9325 break;
9326
9327 case 2:
9331
9335
9339
9343
9345
9348
9351
9354 break;
9355
9356 case 3:
9360
9364
9368
9372
9374
9377
9380
9383 break;
9384 }
9385
9386 // Provide possibility for translation -- may need to add
9387 // face index to this; see ThinLayerBrickOnTetMesh.
9389
9390 // Initialise with identity mapping
9391 for (unsigned i = 0; i < 19; i++)
9392 {
9393 translate[i] = i;
9394 }
9395
9396 // Visit all the nodes on that face
9397 for (unsigned j = 0; j < 19; j++)
9398 {
9399 // Which node is it?
9401
9402 // Get coordinates etc of point from face
9405 Vector<double> x(3);
9408
9409#ifdef PARANOID
9410 // Check that the coordinates match (within tolerance)
9411 double dist = sqrt(pow(brick_node_pt->x(0) - x[0], 2) +
9412 pow(brick_node_pt->x(1) - x[1], 2) +
9413 pow(brick_node_pt->x(2) - x[2], 2));
9415 {
9416 std::ofstream brick0;
9417 std::ofstream brick1;
9418 std::ofstream brick2;
9419 std::ofstream brick3;
9420 brick0.open("full_brick0.dat");
9421 brick1.open("full_brick1.dat");
9422 brick2.open("full_brick2.dat");
9423 brick3.open("full_brick3.dat");
9424 for (unsigned j = 0; j < 27; j++)
9425 {
9426 brick0 << brick_el0_pt->node_pt(j)->x(0) << " "
9427 << brick_el0_pt->node_pt(j)->x(1) << " "
9428 << brick_el0_pt->node_pt(j)->x(2) << "\n";
9429
9430 brick1 << brick_el1_pt->node_pt(j)->x(0) << " "
9431 << brick_el1_pt->node_pt(j)->x(1) << " "
9432 << brick_el1_pt->node_pt(j)->x(2) << "\n";
9433
9434 brick2 << brick_el2_pt->node_pt(j)->x(0) << " "
9435 << brick_el2_pt->node_pt(j)->x(1) << " "
9436 << brick_el2_pt->node_pt(j)->x(2) << "\n";
9437
9438 brick3 << brick_el3_pt->node_pt(j)->x(0) << " "
9439 << brick_el3_pt->node_pt(j)->x(1) << " "
9440 << brick_el3_pt->node_pt(j)->x(2) << "\n";
9441 }
9442 brick0.close();
9443 brick1.close();
9444 brick2.close();
9445 brick3.close();
9446
9447 std::ofstream full_face;
9448 full_face.open("full_face.dat");
9449 for (unsigned j = 0; j < 6; j++)
9450 {
9451 full_face << face_el_pt->node_pt(j)->x(0) << " "
9452 << face_el_pt->node_pt(j)->x(1) << " "
9453 << face_el_pt->node_pt(j)->x(2) << "\n";
9454 }
9455 full_face.close();
9456
9457 // Get normal sign
9458 int normal_sign = face_el_pt->normal_sign();
9459
9460 std::ostringstream error_stream;
9462 << "During assignment of boundary cordinates, the distance\n"
9463 << "between brick node and reference point in \n"
9464 << "triangular FaceElement is " << dist << " which \n"
9465 << "is bigger than the tolerance defined in \n"
9466 << "BrickFromTetMeshHelper::Face_position_tolerance="
9468 << "If this is tolerable, increase the tolerance \n"
9469 << "(it's defined in a namespace and therefore publically\n"
9470 << "accessible). If not, the Face may be inverted in which \n"
9471 << "case you should re-implement the translation scheme,\n"
9472 << "following the pattern used in the "
9473 "ThinLayerBrickOnTetMesh."
9474 << "\nThe required code fragements are already here but \n"
9475 << "the translation is the unit map.\n"
9476 << "To aid the diagnostics, the files full_brick[0-3].dat\n"
9477 << "contain the coordinates of the 27 nodes in the four\n"
9478 << "bricks associated with the current tet and "
9479 "full_face.dat\n"
9480 << "contains the coordinates of the 6 nodes in the "
9481 "FaceElement"
9482 << "\nfrom which the boundary coordinates are extracted.\n"
9483 << "FYI: The normal_sign of the face is: " << normal_sign
9484 << std::endl;
9485 throw OomphLibError(error_stream.str(),
9488 }
9489#endif
9490
9491 // Set boundary stuff
9492 add_boundary_node(b, brick_node_pt);
9493 brick_node_pt->set_coordinates_on_boundary(b, zeta);
9494 }
9495
9496 // Add appropriate brick elements to boundary lookup scheme
9497 switch (face_index)
9498 {
9499 case 0:
9500 Boundary_element_pt[b].push_back(brick_el1_pt);
9501 Face_index_at_boundary[b].push_back(-2);
9502 Boundary_element_pt[b].push_back(brick_el2_pt);
9503 Face_index_at_boundary[b].push_back(-1);
9504 Boundary_element_pt[b].push_back(brick_el3_pt);
9505 Face_index_at_boundary[b].push_back(-2);
9506 break;
9507
9508 case 1:
9509 Boundary_element_pt[b].push_back(brick_el0_pt);
9510 Face_index_at_boundary[b].push_back(-1);
9511 Boundary_element_pt[b].push_back(brick_el2_pt);
9512 Face_index_at_boundary[b].push_back(-2);
9513 Boundary_element_pt[b].push_back(brick_el3_pt);
9514 Face_index_at_boundary[b].push_back(-1);
9515 break;
9516
9517 case 2:
9518 Boundary_element_pt[b].push_back(brick_el0_pt);
9519 Face_index_at_boundary[b].push_back(-3);
9520 Boundary_element_pt[b].push_back(brick_el1_pt);
9521 Face_index_at_boundary[b].push_back(-3);
9522 Boundary_element_pt[b].push_back(brick_el2_pt);
9523 Face_index_at_boundary[b].push_back(-3);
9524 break;
9525
9526 case 3:
9527 Boundary_element_pt[b].push_back(brick_el0_pt);
9528 Face_index_at_boundary[b].push_back(-2);
9529 Boundary_element_pt[b].push_back(brick_el1_pt);
9530 Face_index_at_boundary[b].push_back(-1);
9531 Boundary_element_pt[b].push_back(brick_el3_pt);
9532 Face_index_at_boundary[b].push_back(-3);
9533 break;
9534 }
9535 // Cleanup
9536 delete face_el_pt;
9537 }
9538 }
9539 // Cleanup
9540 delete face_pt;
9541 }
9542 }
9543
9544 // Lookup scheme has now been setup
9545 Lookup_for_elements_next_boundary_is_setup = true;
9546
9547 // Get number of distinct boundaries specified
9548 // in the original xda enumeration.
9549 // unsigned n_xda_boundaries=tet_mesh_pt->nboundary();
9550
9551 // Copy collective IDs across
9552 // Boundary_id.resize(n_xda_boundaries);
9553 // for (unsigned xda_b=0;xda_b<n_xda_boundaries;xda_b++)
9554 // {
9555 // //Boundary_id[xda_b]=tet_mesh_pt->oomph_lib_boundary_ids(xda_b);
9556 // Boundary_id[xda_b]=xda_b;
9557 // }
9558
9559 // Cleanup
9560 for (unsigned e = 0; e < 4; e++)
9561 {
9562 for (unsigned j = 0; j < 8; j++)
9563 {
9564 delete dummy_q_el_pt[e]->node_pt(j);
9565 }
9566 delete dummy_q_el_pt[e];
9567 }
9568 }
9569
9570} // namespace oomph
9571
9572#endif
e
Definition cfortran.h:571
static char t char * s
Definition cfortran.h:568
cstr elem_len * i
Definition cfortran.h:603
void build_mesh(XdaTetMesh< TElement< 3, 3 > > *tet_mesh_pt, TimeStepper *time_stepper_pt)
Build fct: Pass pointer to existing tet mesh.
Dummy QElement to interpolate local coordinates – used in construction of brickified tet mesh.
Definition brick_mesh.h:76
Edge class.
Definition mesh.h:2700
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition elements.h:4342
A general Finite Element class.
Definition elements.h:1317
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
Definition elements.h:1846
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition elements.cc:4705
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition elements.h:2513
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition elements.cc:3992
Node ** Node_pt
Storage for pointers to the nodes in the element.
Definition elements.h:1323
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition elements.h:2542
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
SolidFiniteElement class.
Definition elements.h:3565
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
Triangular Face class.
Definition Telements.h:56
Unstructured tet mesh based on output from Tetgen: http://wias-berlin.de/software/tetgen/.
Definition tetgen_mesh.h:51
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Tet mesh made of quadratic (ten node) tets built from xda input file.
double Face_position_tolerance
Tolerance for mismatch during setup of boundary coordinates.
Definition brick_mesh.cc:37
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).