30 #ifndef OOMPH_BRICK_FROM_TET_MESH_TEMPLATE_CC 31 #define OOMPH_BRICK_FROM_TET_MESH_TEMPLATE_CC 44 template<
class ELEMENT>
47 TimeStepper* time_stepper_pt)
50 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3, 3);
53 bool tet_mesh_is_solid_mesh=
false;
54 if (dynamic_cast<SolidFiniteElement*>(tet_mesh_pt->element_pt(0))!=0)
56 tet_mesh_is_solid_mesh=
true;
63 Vector<Vector<double> > s_face(19);
64 for (
unsigned i=0;i<19;i++)
141 s_face[i][0]=1.0/3.0;
142 s_face[i][1]=1.0/3.0;
149 s_face[i][0]=5.0/24.0;
150 s_face[i][1]=5.0/24.0;
154 s_face[i][0]=5.0/12.0;
155 s_face[i][1]=5.0/12.0;
161 s_face[i][1]=5.0/24.0;
162 s_face[i][0]=7.0/12.0;
166 s_face[i][1]=5.0/12.0;
167 s_face[i][0]=1.0/6.0;
174 s_face[i][0]=5.0/24.0;
175 s_face[i][1]=7.0/12.0;
179 s_face[i][0]=5.0/12.0;
180 s_face[i][1]=1.0/6.0;
187 unsigned nb=tet_mesh_pt->nboundary();
191 Boundary_element_pt.resize(nb);
192 Face_index_at_boundary.resize(nb);
197 std::map<Node*,Node*> tet_node_node_pt;
200 std::map<Edge,Node*> brick_edge_node_pt;
203 std::map<TFace,Node*> tet_face_node_pt;
207 Vector<DummyBrickElement*> dummy_q_el_pt(4);
208 for (
unsigned e=0;e<4;e++)
210 dummy_q_el_pt[e]=
new DummyBrickElement;
211 for (
unsigned j=0;j<8;j++)
213 dummy_q_el_pt[e]->construct_node(j);
218 unsigned n_el_tet=tet_mesh_pt->nelement();
219 for (
unsigned e_tet=0; e_tet<n_el_tet; e_tet++)
222 TElement<3,3>* tet_el_pt=
dynamic_cast<TElement<3,3>*
>(
223 tet_mesh_pt->element_pt(e_tet));
228 std::ostringstream error_stream;
230 <<
"BrickFromTetMesh can only built from tet mesh containing\n" 231 <<
"ten-noded tets.\n";
234 OOMPH_CURRENT_FUNCTION,
235 OOMPH_EXCEPTION_LOCATION);
240 Node* centroid_node_pt=0;
243 Node* top_mid_face_node0_pt=0;
244 Node* right_mid_face_node0_pt=0;
245 Node* back_mid_face_node0_pt=0;
247 Node* top_mid_face_node1_pt=0;
248 Node* right_mid_face_node1_pt=0;
250 Node* top_mid_face_node2_pt=0;
253 FiniteElement* brick_el0_pt=0;
254 FiniteElement* brick_el1_pt=0;
255 FiniteElement* brick_el2_pt=0;
256 FiniteElement* brick_el3_pt=0;
264 for (
unsigned j=0;j<8;j++)
266 Node* nod_pt=dummy_q_el_pt[0]->node_pt(j);
267 Vector<double> s_tet(3);
268 Vector<double> x_tet(3);
272 tet_el_pt->local_coordinate_of_node(0,s_tet);
273 nod_pt->set_value(0,s_tet[0]);
274 nod_pt->set_value(1,s_tet[1]);
275 nod_pt->set_value(2,s_tet[2]);
276 tet_el_pt->interpolated_x(s_tet,x_tet);
277 nod_pt->x(0)=x_tet[0];
278 nod_pt->x(1)=x_tet[1];
279 nod_pt->x(2)=x_tet[2];
282 tet_el_pt->local_coordinate_of_node(4,s_tet);
283 nod_pt->set_value(0,s_tet[0]);
284 nod_pt->set_value(1,s_tet[1]);
285 nod_pt->set_value(2,s_tet[2]);
286 tet_el_pt->interpolated_x(s_tet,x_tet);
287 nod_pt->x(0)=x_tet[0];
288 nod_pt->x(1)=x_tet[1];
289 nod_pt->x(2)=x_tet[2];
292 tet_el_pt->local_coordinate_of_node(6,s_tet);
293 nod_pt->set_value(0,s_tet[0]);
294 nod_pt->set_value(1,s_tet[1]);
295 nod_pt->set_value(2,s_tet[2]);
296 tet_el_pt->interpolated_x(s_tet,x_tet);
297 nod_pt->x(0)=x_tet[0];
298 nod_pt->x(1)=x_tet[1];
299 nod_pt->x(2)=x_tet[2];
307 nod_pt->set_value(0,s_tet[0]);
308 nod_pt->set_value(1,s_tet[1]);
309 nod_pt->set_value(2,s_tet[2]);
310 tet_el_pt->interpolated_x(s_tet,x_tet);
311 nod_pt->x(0)=x_tet[0];
312 nod_pt->x(1)=x_tet[1];
313 nod_pt->x(2)=x_tet[2];
316 tet_el_pt->local_coordinate_of_node(5,s_tet);
317 nod_pt->set_value(0,s_tet[0]);
318 nod_pt->set_value(1,s_tet[1]);
319 nod_pt->set_value(2,s_tet[2]);
320 tet_el_pt->interpolated_x(s_tet,x_tet);
321 nod_pt->x(0)=x_tet[0];
322 nod_pt->x(1)=x_tet[1];
323 nod_pt->x(2)=x_tet[2];
331 nod_pt->set_value(0,s_tet[0]);
332 nod_pt->set_value(1,s_tet[1]);
333 nod_pt->set_value(2,s_tet[2]);
334 tet_el_pt->interpolated_x(s_tet,x_tet);
335 nod_pt->x(0)=x_tet[0];
336 nod_pt->x(1)=x_tet[1];
337 nod_pt->x(2)=x_tet[2];
345 nod_pt->set_value(0,s_tet[0]);
346 nod_pt->set_value(1,s_tet[1]);
347 nod_pt->set_value(2,s_tet[2]);
348 tet_el_pt->interpolated_x(s_tet,x_tet);
349 nod_pt->x(0)=x_tet[0];
350 nod_pt->x(1)=x_tet[1];
351 nod_pt->x(2)=x_tet[2];
358 nod_pt->set_value(0,s_tet[0]);
359 nod_pt->set_value(1,s_tet[1]);
360 nod_pt->set_value(2,s_tet[2]);
361 tet_el_pt->interpolated_x(s_tet,x_tet);
362 nod_pt->x(0)=x_tet[0];
363 nod_pt->x(1)=x_tet[1];
364 nod_pt->x(2)=x_tet[2];
371 FiniteElement* el_pt=
new ELEMENT;
373 Element_pt.push_back(el_pt);
375 TFace face0(tet_el_pt->node_pt(0),
376 tet_el_pt->node_pt(1),
377 tet_el_pt->node_pt(2));
379 TFace face1(tet_el_pt->node_pt(0),
380 tet_el_pt->node_pt(2),
381 tet_el_pt->node_pt(3));
383 TFace face2(tet_el_pt->node_pt(0),
384 tet_el_pt->node_pt(1),
385 tet_el_pt->node_pt(3));
389 Vector<Vector<unsigned> > tet_edge_node(3);
390 tet_edge_node[0].resize(2);
391 tet_edge_node[0][0]=4;
392 tet_edge_node[0][1]=1;
393 tet_edge_node[1].resize(2);
394 tet_edge_node[1][0]=6;
395 tet_edge_node[1][1]=3;
396 tet_edge_node[2].resize(2);
397 tet_edge_node[2][0]=5;
398 tet_edge_node[2][1]=2;
401 unsigned central_tet_vertex=0;
411 tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
412 old_node_pt=tet_node_node_pt[tet_node_pt];
416 if (tet_node_pt->is_on_boundary())
418 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
422 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
424 tet_node_node_pt[tet_node_pt]=new_node_pt;
425 Node_pt.push_back(new_node_pt);
427 Vector<double> s_tet(3);
428 Vector<double> x_tet(3);
429 el_pt->local_coordinate_of_node(j,s);
430 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
431 tet_el_pt->interpolated_x(s_tet,x_tet);
432 new_node_pt->x(0)=x_tet[0];
433 new_node_pt->x(1)=x_tet[1];
434 new_node_pt->x(2)=x_tet[2];
439 el_pt->node_pt(j)=old_node_pt;
449 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
450 old_node_pt=tet_node_node_pt[tet_node_pt];
454 if (tet_node_pt->is_on_boundary())
456 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
460 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
462 tet_node_node_pt[tet_node_pt]=new_node_pt;
463 Node_pt.push_back(new_node_pt);
465 Vector<double> s_tet(3);
466 Vector<double> x_tet(3);
467 el_pt->local_coordinate_of_node(j,s);
468 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
469 tet_el_pt->interpolated_x(s_tet,x_tet);
470 new_node_pt->x(0)=x_tet[0];
471 new_node_pt->x(1)=x_tet[1];
472 new_node_pt->x(2)=x_tet[2];
477 el_pt->node_pt(j)=old_node_pt;
487 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
488 old_node_pt=tet_node_node_pt[tet_node_pt];
492 if (tet_node_pt->is_on_boundary())
494 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
498 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
500 tet_node_node_pt[tet_node_pt]=new_node_pt;
501 Node_pt.push_back(new_node_pt);
503 Vector<double> s_tet(3);
504 Vector<double> x_tet(3);
505 el_pt->local_coordinate_of_node(j,s);
506 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
507 tet_el_pt->interpolated_x(s_tet,x_tet);
508 new_node_pt->x(0)=x_tet[0];
509 new_node_pt->x(1)=x_tet[1];
510 new_node_pt->x(2)=x_tet[2];
515 el_pt->node_pt(j)=old_node_pt;
525 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
526 old_node_pt=tet_node_node_pt[tet_node_pt];
530 if (tet_node_pt->is_on_boundary())
532 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
536 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
538 tet_node_node_pt[tet_node_pt]=new_node_pt;
539 Node_pt.push_back(new_node_pt);
541 Vector<double> s_tet(3);
542 Vector<double> x_tet(3);
543 el_pt->local_coordinate_of_node(j,s);
544 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
545 tet_el_pt->interpolated_x(s_tet,x_tet);
546 new_node_pt->x(0)=x_tet[0];
547 new_node_pt->x(1)=x_tet[1];
548 new_node_pt->x(2)=x_tet[2];
553 el_pt->node_pt(j)=old_node_pt;
565 old_node_pt=tet_face_node_pt[face0];
569 if (face0.is_boundary_face())
571 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
575 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
577 tet_face_node_pt[face0]=new_node_pt;
578 Node_pt.push_back(new_node_pt);
580 Vector<double> s_tet(3);
581 Vector<double> x_tet(3);
582 el_pt->local_coordinate_of_node(j,s);
583 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
584 tet_el_pt->interpolated_x(s_tet,x_tet);
585 new_node_pt->x(0)=x_tet[0];
586 new_node_pt->x(1)=x_tet[1];
587 new_node_pt->x(2)=x_tet[2];
592 el_pt->node_pt(j)=old_node_pt;
602 old_node_pt=tet_face_node_pt[face1];
606 if (face1.is_boundary_face())
608 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
612 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
614 tet_face_node_pt[face1]=new_node_pt;
615 Node_pt.push_back(new_node_pt);
617 Vector<double> s_tet(3);
618 Vector<double> x_tet(3);
619 el_pt->local_coordinate_of_node(j,s);
620 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
621 tet_el_pt->interpolated_x(s_tet,x_tet);
622 new_node_pt->x(0)=x_tet[0];
623 new_node_pt->x(1)=x_tet[1];
624 new_node_pt->x(2)=x_tet[2];
629 el_pt->node_pt(j)=old_node_pt;
639 old_node_pt=tet_face_node_pt[face2];
643 if (face2.is_boundary_face())
645 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
649 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
651 tet_face_node_pt[face2]=new_node_pt;
652 Node_pt.push_back(new_node_pt);
654 Vector<double> s_tet(3);
655 Vector<double> x_tet(3);
656 el_pt->local_coordinate_of_node(j,s);
657 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
658 tet_el_pt->interpolated_x(s_tet,x_tet);
659 new_node_pt->x(0)=x_tet[0];
660 new_node_pt->x(1)=x_tet[1];
661 new_node_pt->x(2)=x_tet[2];
666 el_pt->node_pt(j)=old_node_pt;
677 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
678 centroid_node_pt=new_node_pt;
679 Node_pt.push_back(new_node_pt);
681 Vector<double> s_tet(3);
682 Vector<double> x_tet(3);
683 el_pt->local_coordinate_of_node(j,s);
684 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
685 tet_el_pt->interpolated_x(s_tet,x_tet);
686 new_node_pt->x(0)=x_tet[0];
687 new_node_pt->x(1)=x_tet[1];
688 new_node_pt->x(2)=x_tet[2];
699 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
700 Node_pt.push_back(new_node_pt);
702 Vector<double> s_tet(3);
703 Vector<double> x_tet(3);
704 el_pt->local_coordinate_of_node(j,s);
705 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
706 tet_el_pt->interpolated_x(s_tet,x_tet);
707 new_node_pt->x(0)=x_tet[0];
708 new_node_pt->x(1)=x_tet[1];
709 new_node_pt->x(2)=x_tet[2];
718 Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
719 old_node_pt=brick_edge_node_pt[edge];
723 if (edge.is_boundary_edge())
725 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
729 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
731 brick_edge_node_pt[edge]=new_node_pt;
732 Node_pt.push_back(new_node_pt);
734 Vector<double> s_tet(3);
735 Vector<double> x_tet(3);
736 el_pt->local_coordinate_of_node(j,s);
737 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
738 tet_el_pt->interpolated_x(s_tet,x_tet);
739 new_node_pt->x(0)=x_tet[0];
740 new_node_pt->x(1)=x_tet[1];
741 new_node_pt->x(2)=x_tet[2];
746 el_pt->node_pt(j)=old_node_pt;
756 Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
757 old_node_pt=brick_edge_node_pt[edge];
761 if (edge.is_boundary_edge())
763 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
767 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
769 brick_edge_node_pt[edge]=new_node_pt;
770 Node_pt.push_back(new_node_pt);
772 Vector<double> s_tet(3);
773 Vector<double> x_tet(3);
774 el_pt->local_coordinate_of_node(j,s);
775 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
776 tet_el_pt->interpolated_x(s_tet,x_tet);
777 new_node_pt->x(0)=x_tet[0];
778 new_node_pt->x(1)=x_tet[1];
779 new_node_pt->x(2)=x_tet[2];
784 el_pt->node_pt(j)=old_node_pt;
793 Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
794 old_node_pt=brick_edge_node_pt[edge];
798 if (edge.is_boundary_edge())
800 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
804 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
806 brick_edge_node_pt[edge]=new_node_pt;
807 Node_pt.push_back(new_node_pt);
809 Vector<double> s_tet(3);
810 Vector<double> x_tet(3);
811 el_pt->local_coordinate_of_node(j,s);
812 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
813 tet_el_pt->interpolated_x(s_tet,x_tet);
814 new_node_pt->x(0)=x_tet[0];
815 new_node_pt->x(1)=x_tet[1];
816 new_node_pt->x(2)=x_tet[2];
821 el_pt->node_pt(j)=old_node_pt;
830 Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
831 old_node_pt=brick_edge_node_pt[edge];
835 if (edge.is_boundary_edge())
837 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
841 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
843 brick_edge_node_pt[edge]=new_node_pt;
844 Node_pt.push_back(new_node_pt);
846 Vector<double> s_tet(3);
847 Vector<double> x_tet(3);
848 el_pt->local_coordinate_of_node(j,s);
849 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
850 tet_el_pt->interpolated_x(s_tet,x_tet);
851 new_node_pt->x(0)=x_tet[0];
852 new_node_pt->x(1)=x_tet[1];
853 new_node_pt->x(2)=x_tet[2];
858 el_pt->node_pt(j)=old_node_pt;
867 Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
868 old_node_pt=brick_edge_node_pt[edge];
872 if (edge.is_boundary_edge())
874 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
878 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
880 brick_edge_node_pt[edge]=new_node_pt;
881 Node_pt.push_back(new_node_pt);
883 Vector<double> s_tet(3);
884 Vector<double> x_tet(3);
885 el_pt->local_coordinate_of_node(j,s);
886 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
887 tet_el_pt->interpolated_x(s_tet,x_tet);
888 new_node_pt->x(0)=x_tet[0];
889 new_node_pt->x(1)=x_tet[1];
890 new_node_pt->x(2)=x_tet[2];
895 el_pt->node_pt(j)=old_node_pt;
905 Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
906 old_node_pt=brick_edge_node_pt[edge];
910 if (edge.is_boundary_edge())
912 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
916 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
918 brick_edge_node_pt[edge]=new_node_pt;
919 Node_pt.push_back(new_node_pt);
921 Vector<double> s_tet(3);
922 Vector<double> x_tet(3);
923 el_pt->local_coordinate_of_node(j,s);
924 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
925 tet_el_pt->interpolated_x(s_tet,x_tet);
926 new_node_pt->x(0)=x_tet[0];
927 new_node_pt->x(1)=x_tet[1];
928 new_node_pt->x(2)=x_tet[2];
933 el_pt->node_pt(j)=old_node_pt;
942 Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
943 old_node_pt=brick_edge_node_pt[edge];
947 if (edge.is_boundary_edge())
949 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
953 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
955 brick_edge_node_pt[edge]=new_node_pt;
956 Node_pt.push_back(new_node_pt);
958 Vector<double> s_tet(3);
959 Vector<double> x_tet(3);
960 el_pt->local_coordinate_of_node(j,s);
961 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
962 tet_el_pt->interpolated_x(s_tet,x_tet);
963 new_node_pt->x(0)=x_tet[0];
964 new_node_pt->x(1)=x_tet[1];
965 new_node_pt->x(2)=x_tet[2];
970 el_pt->node_pt(j)=old_node_pt;
980 Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
981 old_node_pt=brick_edge_node_pt[edge];
985 if (edge.is_boundary_edge())
987 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
991 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
993 brick_edge_node_pt[edge]=new_node_pt;
994 Node_pt.push_back(new_node_pt);
996 Vector<double> s_tet(3);
997 Vector<double> x_tet(3);
998 el_pt->local_coordinate_of_node(j,s);
999 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1000 tet_el_pt->interpolated_x(s_tet,x_tet);
1001 new_node_pt->x(0)=x_tet[0];
1002 new_node_pt->x(1)=x_tet[1];
1003 new_node_pt->x(2)=x_tet[2];
1008 el_pt->node_pt(j)=old_node_pt;
1017 Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
1018 old_node_pt=brick_edge_node_pt[edge];
1021 Node* new_node_pt=0;
1022 if (edge.is_boundary_edge())
1024 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1028 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1030 brick_edge_node_pt[edge]=new_node_pt;
1031 Node_pt.push_back(new_node_pt);
1032 Vector<double> s(3);
1033 Vector<double> s_tet(3);
1034 Vector<double> x_tet(3);
1035 el_pt->local_coordinate_of_node(j,s);
1036 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1037 tet_el_pt->interpolated_x(s_tet,x_tet);
1038 new_node_pt->x(0)=x_tet[0];
1039 new_node_pt->x(1)=x_tet[1];
1040 new_node_pt->x(2)=x_tet[2];
1045 el_pt->node_pt(j)=old_node_pt;
1055 Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
1056 old_node_pt=brick_edge_node_pt[edge];
1059 Node* new_node_pt=0;
1060 if (edge.is_boundary_edge())
1062 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1066 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1068 brick_edge_node_pt[edge]=new_node_pt;
1069 Node_pt.push_back(new_node_pt);
1070 Vector<double> s(3);
1071 Vector<double> s_tet(3);
1072 Vector<double> x_tet(3);
1073 el_pt->local_coordinate_of_node(j,s);
1074 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1075 tet_el_pt->interpolated_x(s_tet,x_tet);
1076 new_node_pt->x(0)=x_tet[0];
1077 new_node_pt->x(1)=x_tet[1];
1078 new_node_pt->x(2)=x_tet[2];
1083 el_pt->node_pt(j)=old_node_pt;
1093 Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
1094 old_node_pt=brick_edge_node_pt[edge];
1097 Node* new_node_pt=0;
1098 if (edge.is_boundary_edge())
1100 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1104 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1106 brick_edge_node_pt[edge]=new_node_pt;
1107 Node_pt.push_back(new_node_pt);
1108 Vector<double> s(3);
1109 Vector<double> s_tet(3);
1110 Vector<double> x_tet(3);
1111 el_pt->local_coordinate_of_node(j,s);
1112 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1113 tet_el_pt->interpolated_x(s_tet,x_tet);
1114 new_node_pt->x(0)=x_tet[0];
1115 new_node_pt->x(1)=x_tet[1];
1116 new_node_pt->x(2)=x_tet[2];
1121 el_pt->node_pt(j)=old_node_pt;
1131 Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
1132 old_node_pt=brick_edge_node_pt[edge];
1135 Node* new_node_pt=0;
1136 if (edge.is_boundary_edge())
1138 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1142 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1144 brick_edge_node_pt[edge]=new_node_pt;
1145 Node_pt.push_back(new_node_pt);
1146 Vector<double> s(3);
1147 Vector<double> s_tet(3);
1148 Vector<double> x_tet(3);
1149 el_pt->local_coordinate_of_node(j,s);
1150 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1151 tet_el_pt->interpolated_x(s_tet,x_tet);
1152 new_node_pt->x(0)=x_tet[0];
1153 new_node_pt->x(1)=x_tet[1];
1154 new_node_pt->x(2)=x_tet[2];
1159 el_pt->node_pt(j)=old_node_pt;
1170 TFace face(tet_el_pt->node_pt( central_tet_vertex),
1171 tet_el_pt->node_pt(tet_edge_node[0][0]),
1172 tet_el_pt->node_pt(tet_edge_node[2][0]));
1174 old_node_pt=tet_face_node_pt[face];
1177 Node* new_node_pt=0;
1178 if (face.is_boundary_face())
1180 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1184 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1186 tet_face_node_pt[face]=new_node_pt;
1187 Node_pt.push_back(new_node_pt);
1188 Vector<double> s(3);
1189 Vector<double> s_tet(3);
1190 Vector<double> x_tet(3);
1191 el_pt->local_coordinate_of_node(j,s);
1192 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1193 tet_el_pt->interpolated_x(s_tet,x_tet);
1194 new_node_pt->x(0)=x_tet[0];
1195 new_node_pt->x(1)=x_tet[1];
1196 new_node_pt->x(2)=x_tet[2];
1201 el_pt->node_pt(j)=old_node_pt;
1213 TFace face(tet_el_pt->node_pt(central_tet_vertex),
1214 tet_el_pt->node_pt(tet_edge_node[1][0]),
1215 tet_el_pt->node_pt(tet_edge_node[2][0]));
1217 old_node_pt=tet_face_node_pt[face];
1220 Node* new_node_pt=0;
1221 if (face.is_boundary_face())
1223 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1227 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1229 tet_face_node_pt[face]=new_node_pt;
1230 Node_pt.push_back(new_node_pt);
1231 Vector<double> s(3);
1232 Vector<double> s_tet(3);
1233 Vector<double> x_tet(3);
1234 el_pt->local_coordinate_of_node(j,s);
1235 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1236 tet_el_pt->interpolated_x(s_tet,x_tet);
1237 new_node_pt->x(0)=x_tet[0];
1238 new_node_pt->x(1)=x_tet[1];
1239 new_node_pt->x(2)=x_tet[2];
1244 el_pt->node_pt(j)=old_node_pt;
1256 TFace face(tet_el_pt->node_pt(central_tet_vertex),
1257 tet_el_pt->node_pt(tet_edge_node[0][0]),
1258 tet_el_pt->node_pt(tet_edge_node[1][0]));
1260 old_node_pt=tet_face_node_pt[face];
1263 Node* new_node_pt=0;
1264 if (face.is_boundary_face())
1266 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1270 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1272 tet_face_node_pt[face]=new_node_pt;
1273 Node_pt.push_back(new_node_pt);
1274 Vector<double> s(3);
1275 Vector<double> s_tet(3);
1276 Vector<double> x_tet(3);
1277 el_pt->local_coordinate_of_node(j,s);
1278 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1279 tet_el_pt->interpolated_x(s_tet,x_tet);
1280 new_node_pt->x(0)=x_tet[0];
1281 new_node_pt->x(1)=x_tet[1];
1282 new_node_pt->x(2)=x_tet[2];
1287 el_pt->node_pt(j)=old_node_pt;
1295 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1296 Node_pt.push_back(new_node_pt);
1297 Vector<double> s(3);
1298 Vector<double> s_tet(3);
1299 Vector<double> x_tet(3);
1300 el_pt->local_coordinate_of_node(j,s);
1301 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1302 top_mid_face_node0_pt=new_node_pt;
1303 tet_el_pt->interpolated_x(s_tet,x_tet);
1304 new_node_pt->x(0)=x_tet[0];
1305 new_node_pt->x(1)=x_tet[1];
1306 new_node_pt->x(2)=x_tet[2];
1314 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1315 Node_pt.push_back(new_node_pt);
1316 Vector<double> s(3);
1317 Vector<double> s_tet(3);
1318 Vector<double> x_tet(3);
1319 el_pt->local_coordinate_of_node(j,s);
1320 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1321 right_mid_face_node0_pt=new_node_pt;
1322 tet_el_pt->interpolated_x(s_tet,x_tet);
1323 new_node_pt->x(0)=x_tet[0];
1324 new_node_pt->x(1)=x_tet[1];
1325 new_node_pt->x(2)=x_tet[2];
1332 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1333 Node_pt.push_back(new_node_pt);
1334 Vector<double> s(3);
1335 Vector<double> s_tet(3);
1336 Vector<double> x_tet(3);
1337 el_pt->local_coordinate_of_node(j,s);
1338 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
1339 back_mid_face_node0_pt=new_node_pt;
1340 tet_el_pt->interpolated_x(s_tet,x_tet);
1341 new_node_pt->x(0)=x_tet[0];
1342 new_node_pt->x(1)=x_tet[1];
1343 new_node_pt->x(2)=x_tet[2];
1353 for (
unsigned j=0;j<8;j++)
1355 Node* nod_pt=dummy_q_el_pt[1]->node_pt(j);
1356 Vector<double> s_tet(3);
1357 Vector<double> x_tet(3);
1361 tet_el_pt->local_coordinate_of_node(1,s_tet);
1362 nod_pt->set_value(0,s_tet[0]);
1363 nod_pt->set_value(1,s_tet[1]);
1364 nod_pt->set_value(2,s_tet[2]);
1365 tet_el_pt->interpolated_x(s_tet,x_tet);
1366 nod_pt->x(0)=x_tet[0];
1367 nod_pt->x(1)=x_tet[1];
1368 nod_pt->x(2)=x_tet[2];
1371 tet_el_pt->local_coordinate_of_node(9,s_tet);
1372 nod_pt->set_value(0,s_tet[0]);
1373 nod_pt->set_value(1,s_tet[1]);
1374 nod_pt->set_value(2,s_tet[2]);
1375 tet_el_pt->interpolated_x(s_tet,x_tet);
1376 nod_pt->x(0)=x_tet[0];
1377 nod_pt->x(1)=x_tet[1];
1378 nod_pt->x(2)=x_tet[2];
1381 tet_el_pt->local_coordinate_of_node(4,s_tet);
1382 nod_pt->set_value(0,s_tet[0]);
1383 nod_pt->set_value(1,s_tet[1]);
1384 nod_pt->set_value(2,s_tet[2]);
1385 tet_el_pt->interpolated_x(s_tet,x_tet);
1386 nod_pt->x(0)=x_tet[0];
1387 nod_pt->x(1)=x_tet[1];
1388 nod_pt->x(2)=x_tet[2];
1396 nod_pt->set_value(0,s_tet[0]);
1397 nod_pt->set_value(1,s_tet[1]);
1398 nod_pt->set_value(2,s_tet[2]);
1399 tet_el_pt->interpolated_x(s_tet,x_tet);
1400 nod_pt->x(0)=x_tet[0];
1401 nod_pt->x(1)=x_tet[1];
1402 nod_pt->x(2)=x_tet[2];
1405 tet_el_pt->local_coordinate_of_node(7,s_tet);
1406 nod_pt->set_value(0,s_tet[0]);
1407 nod_pt->set_value(1,s_tet[1]);
1408 nod_pt->set_value(2,s_tet[2]);
1409 tet_el_pt->interpolated_x(s_tet,x_tet);
1410 nod_pt->x(0)=x_tet[0];
1411 nod_pt->x(1)=x_tet[1];
1412 nod_pt->x(2)=x_tet[2];
1420 nod_pt->set_value(0,s_tet[0]);
1421 nod_pt->set_value(1,s_tet[1]);
1422 nod_pt->set_value(2,s_tet[2]);
1423 tet_el_pt->interpolated_x(s_tet,x_tet);
1424 nod_pt->x(0)=x_tet[0];
1425 nod_pt->x(1)=x_tet[1];
1426 nod_pt->x(2)=x_tet[2];
1434 nod_pt->set_value(0,s_tet[0]);
1435 nod_pt->set_value(1,s_tet[1]);
1436 nod_pt->set_value(2,s_tet[2]);
1437 tet_el_pt->interpolated_x(s_tet,x_tet);
1438 nod_pt->x(0)=x_tet[0];
1439 nod_pt->x(1)=x_tet[1];
1440 nod_pt->x(2)=x_tet[2];
1447 nod_pt->set_value(0,s_tet[0]);
1448 nod_pt->set_value(1,s_tet[1]);
1449 nod_pt->set_value(2,s_tet[2]);
1450 tet_el_pt->interpolated_x(s_tet,x_tet);
1451 nod_pt->x(0)=x_tet[0];
1452 nod_pt->x(1)=x_tet[1];
1453 nod_pt->x(2)=x_tet[2];
1460 FiniteElement* el_pt=
new ELEMENT;
1462 Element_pt.push_back(el_pt);
1464 TFace face0(tet_el_pt->node_pt(1),
1465 tet_el_pt->node_pt(3),
1466 tet_el_pt->node_pt(2));
1468 TFace face1(tet_el_pt->node_pt(1),
1469 tet_el_pt->node_pt(0),
1470 tet_el_pt->node_pt(2));
1472 TFace face2(tet_el_pt->node_pt(1),
1473 tet_el_pt->node_pt(0),
1474 tet_el_pt->node_pt(3));
1477 Vector<Vector<unsigned> > tet_edge_node(3);
1478 tet_edge_node[0].resize(2);
1479 tet_edge_node[0][0]=9;
1480 tet_edge_node[0][1]=3;
1481 tet_edge_node[1].resize(2);
1482 tet_edge_node[1][0]=4;
1483 tet_edge_node[1][1]=0;
1484 tet_edge_node[2].resize(2);
1485 tet_edge_node[2][0]=7;
1486 tet_edge_node[2][1]=2;
1489 unsigned central_tet_vertex=1;
1491 Node* tet_node_pt=0;
1492 Node* old_node_pt=0;
1499 tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
1500 old_node_pt=tet_node_node_pt[tet_node_pt];
1503 Node* new_node_pt=0;
1504 if (tet_node_pt->is_on_boundary())
1506 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1510 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1512 tet_node_node_pt[tet_node_pt]=new_node_pt;
1513 Node_pt.push_back(new_node_pt);
1514 Vector<double> s(3);
1515 Vector<double> s_tet(3);
1516 Vector<double> x_tet(3);
1517 el_pt->local_coordinate_of_node(j,s);
1518 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1519 tet_el_pt->interpolated_x(s_tet,x_tet);
1520 new_node_pt->x(0)=x_tet[0];
1521 new_node_pt->x(1)=x_tet[1];
1522 new_node_pt->x(2)=x_tet[2];
1527 el_pt->node_pt(j)=old_node_pt;
1537 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
1538 old_node_pt=tet_node_node_pt[tet_node_pt];
1541 Node* new_node_pt=0;
1542 if (tet_node_pt->is_on_boundary())
1544 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1548 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1550 tet_node_node_pt[tet_node_pt]=new_node_pt;
1551 Node_pt.push_back(new_node_pt);
1552 Vector<double> s(3);
1553 Vector<double> s_tet(3);
1554 Vector<double> x_tet(3);
1555 el_pt->local_coordinate_of_node(j,s);
1556 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1557 tet_el_pt->interpolated_x(s_tet,x_tet);
1558 new_node_pt->x(0)=x_tet[0];
1559 new_node_pt->x(1)=x_tet[1];
1560 new_node_pt->x(2)=x_tet[2];
1565 el_pt->node_pt(j)=old_node_pt;
1575 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
1576 old_node_pt=tet_node_node_pt[tet_node_pt];
1579 Node* new_node_pt=0;
1580 if (tet_node_pt->is_on_boundary())
1582 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1586 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1588 tet_node_node_pt[tet_node_pt]=new_node_pt;
1589 Node_pt.push_back(new_node_pt);
1590 Vector<double> s(3);
1591 Vector<double> s_tet(3);
1592 Vector<double> x_tet(3);
1593 el_pt->local_coordinate_of_node(j,s);
1594 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1595 tet_el_pt->interpolated_x(s_tet,x_tet);
1596 new_node_pt->x(0)=x_tet[0];
1597 new_node_pt->x(1)=x_tet[1];
1598 new_node_pt->x(2)=x_tet[2];
1603 el_pt->node_pt(j)=old_node_pt;
1613 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
1614 old_node_pt=tet_node_node_pt[tet_node_pt];
1617 Node* new_node_pt=0;
1618 if (tet_node_pt->is_on_boundary())
1620 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1624 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1626 tet_node_node_pt[tet_node_pt]=new_node_pt;
1627 Node_pt.push_back(new_node_pt);
1628 Vector<double> s(3);
1629 Vector<double> s_tet(3);
1630 Vector<double> x_tet(3);
1631 el_pt->local_coordinate_of_node(j,s);
1632 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1633 tet_el_pt->interpolated_x(s_tet,x_tet);
1634 new_node_pt->x(0)=x_tet[0];
1635 new_node_pt->x(1)=x_tet[1];
1636 new_node_pt->x(2)=x_tet[2];
1641 el_pt->node_pt(j)=old_node_pt;
1652 old_node_pt=tet_face_node_pt[face0];
1655 Node* new_node_pt=0;
1656 if (face0.is_boundary_face())
1658 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1662 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1664 tet_face_node_pt[face0]=new_node_pt;
1665 Node_pt.push_back(new_node_pt);
1666 Vector<double> s(3);
1667 Vector<double> s_tet(3);
1668 Vector<double> x_tet(3);
1669 el_pt->local_coordinate_of_node(j,s);
1670 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1671 tet_el_pt->interpolated_x(s_tet,x_tet);
1672 new_node_pt->x(0)=x_tet[0];
1673 new_node_pt->x(1)=x_tet[1];
1674 new_node_pt->x(2)=x_tet[2];
1679 el_pt->node_pt(j)=old_node_pt;
1688 old_node_pt=tet_face_node_pt[face1];
1691 Node* new_node_pt=0;
1692 if (face1.is_boundary_face())
1694 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1698 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1700 tet_face_node_pt[face1]=new_node_pt;
1701 Node_pt.push_back(new_node_pt);
1702 Vector<double> s(3);
1703 Vector<double> s_tet(3);
1704 Vector<double> x_tet(3);
1705 el_pt->local_coordinate_of_node(j,s);
1706 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1707 tet_el_pt->interpolated_x(s_tet,x_tet);
1708 new_node_pt->x(0)=x_tet[0];
1709 new_node_pt->x(1)=x_tet[1];
1710 new_node_pt->x(2)=x_tet[2];
1715 el_pt->node_pt(j)=old_node_pt;
1724 old_node_pt=tet_face_node_pt[face2];
1727 Node* new_node_pt=0;
1728 if (face2.is_boundary_face())
1730 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1734 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1736 tet_face_node_pt[face2]=new_node_pt;
1737 Node_pt.push_back(new_node_pt);
1738 Vector<double> s(3);
1739 Vector<double> s_tet(3);
1740 Vector<double> x_tet(3);
1741 el_pt->local_coordinate_of_node(j,s);
1742 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1743 tet_el_pt->interpolated_x(s_tet,x_tet);
1744 new_node_pt->x(0)=x_tet[0];
1745 new_node_pt->x(1)=x_tet[1];
1746 new_node_pt->x(2)=x_tet[2];
1751 el_pt->node_pt(j)=old_node_pt;
1762 el_pt->node_pt(j)=centroid_node_pt;
1772 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1773 Node_pt.push_back(new_node_pt);
1774 Vector<double> s(3);
1775 Vector<double> s_tet(3);
1776 Vector<double> x_tet(3);
1777 el_pt->local_coordinate_of_node(j,s);
1778 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1779 tet_el_pt->interpolated_x(s_tet,x_tet);
1780 new_node_pt->x(0)=x_tet[0];
1781 new_node_pt->x(1)=x_tet[1];
1782 new_node_pt->x(2)=x_tet[2];
1792 Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
1793 old_node_pt=brick_edge_node_pt[edge];
1796 Node* new_node_pt=0;
1797 if (edge.is_boundary_edge())
1799 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1803 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1805 brick_edge_node_pt[edge]=new_node_pt;
1806 Node_pt.push_back(new_node_pt);
1807 Vector<double> s(3);
1808 Vector<double> s_tet(3);
1809 Vector<double> x_tet(3);
1810 el_pt->local_coordinate_of_node(j,s);
1811 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1812 tet_el_pt->interpolated_x(s_tet,x_tet);
1813 new_node_pt->x(0)=x_tet[0];
1814 new_node_pt->x(1)=x_tet[1];
1815 new_node_pt->x(2)=x_tet[2];
1820 el_pt->node_pt(j)=old_node_pt;
1830 Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
1831 old_node_pt=brick_edge_node_pt[edge];
1834 Node* new_node_pt=0;
1835 if (edge.is_boundary_edge())
1837 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1841 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1843 brick_edge_node_pt[edge]=new_node_pt;
1844 Node_pt.push_back(new_node_pt);
1845 Vector<double> s(3);
1846 Vector<double> s_tet(3);
1847 Vector<double> x_tet(3);
1848 el_pt->local_coordinate_of_node(j,s);
1849 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1850 tet_el_pt->interpolated_x(s_tet,x_tet);
1851 new_node_pt->x(0)=x_tet[0];
1852 new_node_pt->x(1)=x_tet[1];
1853 new_node_pt->x(2)=x_tet[2];
1858 el_pt->node_pt(j)=old_node_pt;
1868 Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
1869 old_node_pt=brick_edge_node_pt[edge];
1872 Node* new_node_pt=0;
1873 if (edge.is_boundary_edge())
1875 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1879 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1881 brick_edge_node_pt[edge]=new_node_pt;
1882 Node_pt.push_back(new_node_pt);
1883 Vector<double> s(3);
1884 Vector<double> s_tet(3);
1885 Vector<double> x_tet(3);
1886 el_pt->local_coordinate_of_node(j,s);
1887 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1888 tet_el_pt->interpolated_x(s_tet,x_tet);
1889 new_node_pt->x(0)=x_tet[0];
1890 new_node_pt->x(1)=x_tet[1];
1891 new_node_pt->x(2)=x_tet[2];
1896 el_pt->node_pt(j)=old_node_pt;
1905 Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
1906 old_node_pt=brick_edge_node_pt[edge];
1909 Node* new_node_pt=0;
1910 if (edge.is_boundary_edge())
1912 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1916 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1918 brick_edge_node_pt[edge]=new_node_pt;
1919 Node_pt.push_back(new_node_pt);
1920 Vector<double> s(3);
1921 Vector<double> s_tet(3);
1922 Vector<double> x_tet(3);
1923 el_pt->local_coordinate_of_node(j,s);
1924 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1925 tet_el_pt->interpolated_x(s_tet,x_tet);
1926 new_node_pt->x(0)=x_tet[0];
1927 new_node_pt->x(1)=x_tet[1];
1928 new_node_pt->x(2)=x_tet[2];
1933 el_pt->node_pt(j)=old_node_pt;
1942 Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
1943 old_node_pt=brick_edge_node_pt[edge];
1946 Node* new_node_pt=0;
1947 if (edge.is_boundary_edge())
1949 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1953 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1955 brick_edge_node_pt[edge]=new_node_pt;
1956 Node_pt.push_back(new_node_pt);
1957 Vector<double> s(3);
1958 Vector<double> s_tet(3);
1959 Vector<double> x_tet(3);
1960 el_pt->local_coordinate_of_node(j,s);
1961 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
1962 tet_el_pt->interpolated_x(s_tet,x_tet);
1963 new_node_pt->x(0)=x_tet[0];
1964 new_node_pt->x(1)=x_tet[1];
1965 new_node_pt->x(2)=x_tet[2];
1970 el_pt->node_pt(j)=old_node_pt;
1980 Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
1981 old_node_pt=brick_edge_node_pt[edge];
1984 Node* new_node_pt=0;
1985 if (edge.is_boundary_edge())
1987 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
1991 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
1993 brick_edge_node_pt[edge]=new_node_pt;
1994 Node_pt.push_back(new_node_pt);
1995 Vector<double> s(3);
1996 Vector<double> s_tet(3);
1997 Vector<double> x_tet(3);
1998 el_pt->local_coordinate_of_node(j,s);
1999 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2000 tet_el_pt->interpolated_x(s_tet,x_tet);
2001 new_node_pt->x(0)=x_tet[0];
2002 new_node_pt->x(1)=x_tet[1];
2003 new_node_pt->x(2)=x_tet[2];
2008 el_pt->node_pt(j)=old_node_pt;
2017 Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
2018 old_node_pt=brick_edge_node_pt[edge];
2021 Node* new_node_pt=0;
2022 if (edge.is_boundary_edge())
2024 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2028 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2030 brick_edge_node_pt[edge]=new_node_pt;
2031 Node_pt.push_back(new_node_pt);
2032 Vector<double> s(3);
2033 Vector<double> s_tet(3);
2034 Vector<double> x_tet(3);
2035 el_pt->local_coordinate_of_node(j,s);
2036 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2037 tet_el_pt->interpolated_x(s_tet,x_tet);
2038 new_node_pt->x(0)=x_tet[0];
2039 new_node_pt->x(1)=x_tet[1];
2040 new_node_pt->x(2)=x_tet[2];
2045 el_pt->node_pt(j)=old_node_pt;
2055 Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
2056 old_node_pt=brick_edge_node_pt[edge];
2059 Node* new_node_pt=0;
2060 if (edge.is_boundary_edge())
2062 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2066 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2068 brick_edge_node_pt[edge]=new_node_pt;
2069 Node_pt.push_back(new_node_pt);
2070 Vector<double> s(3);
2071 Vector<double> s_tet(3);
2072 Vector<double> x_tet(3);
2073 el_pt->local_coordinate_of_node(j,s);
2074 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2075 tet_el_pt->interpolated_x(s_tet,x_tet);
2076 new_node_pt->x(0)=x_tet[0];
2077 new_node_pt->x(1)=x_tet[1];
2078 new_node_pt->x(2)=x_tet[2];
2083 el_pt->node_pt(j)=old_node_pt;
2092 Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
2093 old_node_pt=brick_edge_node_pt[edge];
2096 Node* new_node_pt=0;
2097 if (edge.is_boundary_edge())
2099 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2103 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2105 brick_edge_node_pt[edge]=new_node_pt;
2106 Node_pt.push_back(new_node_pt);
2107 Vector<double> s(3);
2108 Vector<double> s_tet(3);
2109 Vector<double> x_tet(3);
2110 el_pt->local_coordinate_of_node(j,s);
2111 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2112 tet_el_pt->interpolated_x(s_tet,x_tet);
2113 new_node_pt->x(0)=x_tet[0];
2114 new_node_pt->x(1)=x_tet[1];
2115 new_node_pt->x(2)=x_tet[2];
2120 el_pt->node_pt(j)=old_node_pt;
2130 Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
2131 old_node_pt=brick_edge_node_pt[edge];
2134 Node* new_node_pt=0;
2135 if (edge.is_boundary_edge())
2137 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2141 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2143 brick_edge_node_pt[edge]=new_node_pt;
2144 Node_pt.push_back(new_node_pt);
2145 Vector<double> s(3);
2146 Vector<double> s_tet(3);
2147 Vector<double> x_tet(3);
2148 el_pt->local_coordinate_of_node(j,s);
2149 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2150 tet_el_pt->interpolated_x(s_tet,x_tet);
2151 new_node_pt->x(0)=x_tet[0];
2152 new_node_pt->x(1)=x_tet[1];
2153 new_node_pt->x(2)=x_tet[2];
2158 el_pt->node_pt(j)=old_node_pt;
2168 Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
2169 old_node_pt=brick_edge_node_pt[edge];
2172 Node* new_node_pt=0;
2173 if (edge.is_boundary_edge())
2175 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2179 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2181 brick_edge_node_pt[edge]=new_node_pt;
2182 Node_pt.push_back(new_node_pt);
2183 Vector<double> s(3);
2184 Vector<double> s_tet(3);
2185 Vector<double> x_tet(3);
2186 el_pt->local_coordinate_of_node(j,s);
2187 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2188 tet_el_pt->interpolated_x(s_tet,x_tet);
2189 new_node_pt->x(0)=x_tet[0];
2190 new_node_pt->x(1)=x_tet[1];
2191 new_node_pt->x(2)=x_tet[2];
2196 el_pt->node_pt(j)=old_node_pt;
2206 Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
2207 old_node_pt=brick_edge_node_pt[edge];
2210 Node* new_node_pt=0;
2211 if (edge.is_boundary_edge())
2213 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2217 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2219 brick_edge_node_pt[edge]=new_node_pt;
2220 Node_pt.push_back(new_node_pt);
2221 Vector<double> s(3);
2222 Vector<double> s_tet(3);
2223 Vector<double> x_tet(3);
2224 el_pt->local_coordinate_of_node(j,s);
2225 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2226 tet_el_pt->interpolated_x(s_tet,x_tet);
2227 new_node_pt->x(0)=x_tet[0];
2228 new_node_pt->x(1)=x_tet[1];
2229 new_node_pt->x(2)=x_tet[2];
2234 el_pt->node_pt(j)=old_node_pt;
2245 TFace face(tet_el_pt->node_pt( central_tet_vertex),
2246 tet_el_pt->node_pt(tet_edge_node[0][0]),
2247 tet_el_pt->node_pt(tet_edge_node[2][0]));
2249 old_node_pt=tet_face_node_pt[face];
2252 Node* new_node_pt=0;
2253 if (face.is_boundary_face())
2255 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2259 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2261 tet_face_node_pt[face]=new_node_pt;
2262 Node_pt.push_back(new_node_pt);
2263 Vector<double> s(3);
2264 Vector<double> s_tet(3);
2265 Vector<double> x_tet(3);
2266 el_pt->local_coordinate_of_node(j,s);
2267 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2268 tet_el_pt->interpolated_x(s_tet,x_tet);
2269 new_node_pt->x(0)=x_tet[0];
2270 new_node_pt->x(1)=x_tet[1];
2271 new_node_pt->x(2)=x_tet[2];
2276 el_pt->node_pt(j)=old_node_pt;
2288 TFace face(tet_el_pt->node_pt(central_tet_vertex),
2289 tet_el_pt->node_pt(tet_edge_node[1][0]),
2290 tet_el_pt->node_pt(tet_edge_node[2][0]));
2292 old_node_pt=tet_face_node_pt[face];
2295 Node* new_node_pt=0;
2296 if (face.is_boundary_face())
2298 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2302 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2304 tet_face_node_pt[face]=new_node_pt;
2305 Node_pt.push_back(new_node_pt);
2306 Vector<double> s(3);
2307 Vector<double> s_tet(3);
2308 Vector<double> x_tet(3);
2309 el_pt->local_coordinate_of_node(j,s);
2310 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2311 tet_el_pt->interpolated_x(s_tet,x_tet);
2312 new_node_pt->x(0)=x_tet[0];
2313 new_node_pt->x(1)=x_tet[1];
2314 new_node_pt->x(2)=x_tet[2];
2319 el_pt->node_pt(j)=old_node_pt;
2331 TFace face(tet_el_pt->node_pt(central_tet_vertex),
2332 tet_el_pt->node_pt(tet_edge_node[0][0]),
2333 tet_el_pt->node_pt(tet_edge_node[1][0]));
2335 old_node_pt=tet_face_node_pt[face];
2338 Node* new_node_pt=0;
2339 if (face.is_boundary_face())
2341 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2345 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2347 tet_face_node_pt[face]=new_node_pt;
2348 Node_pt.push_back(new_node_pt);
2349 Vector<double> s(3);
2350 Vector<double> s_tet(3);
2351 Vector<double> x_tet(3);
2352 el_pt->local_coordinate_of_node(j,s);
2353 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2354 tet_el_pt->interpolated_x(s_tet,x_tet);
2355 new_node_pt->x(0)=x_tet[0];
2356 new_node_pt->x(1)=x_tet[1];
2357 new_node_pt->x(2)=x_tet[2];
2362 el_pt->node_pt(j)=old_node_pt;
2370 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2371 Node_pt.push_back(new_node_pt);
2372 Vector<double> s(3);
2373 Vector<double> s_tet(3);
2374 Vector<double> x_tet(3);
2375 el_pt->local_coordinate_of_node(j,s);
2376 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2377 top_mid_face_node1_pt=new_node_pt;
2378 tet_el_pt->interpolated_x(s_tet,x_tet);
2379 new_node_pt->x(0)=x_tet[0];
2380 new_node_pt->x(1)=x_tet[1];
2381 new_node_pt->x(2)=x_tet[2];
2389 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2390 Node_pt.push_back(new_node_pt);
2391 Vector<double> s(3);
2392 Vector<double> s_tet(3);
2393 Vector<double> x_tet(3);
2394 el_pt->local_coordinate_of_node(j,s);
2395 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
2396 right_mid_face_node1_pt=new_node_pt;
2397 tet_el_pt->interpolated_x(s_tet,x_tet);
2398 new_node_pt->x(0)=x_tet[0];
2399 new_node_pt->x(1)=x_tet[1];
2400 new_node_pt->x(2)=x_tet[2];
2409 el_pt->node_pt(j)=right_mid_face_node0_pt;
2418 for (
unsigned j=0;j<8;j++)
2420 Node* nod_pt=dummy_q_el_pt[2]->node_pt(j);
2421 Vector<double> s_tet(3);
2422 Vector<double> x_tet(3);
2426 tet_el_pt->local_coordinate_of_node(3,s_tet);
2427 nod_pt->set_value(0,s_tet[0]);
2428 nod_pt->set_value(1,s_tet[1]);
2429 nod_pt->set_value(2,s_tet[2]);
2430 tet_el_pt->interpolated_x(s_tet,x_tet);
2431 nod_pt->x(0)=x_tet[0];
2432 nod_pt->x(1)=x_tet[1];
2433 nod_pt->x(2)=x_tet[2];
2436 tet_el_pt->local_coordinate_of_node(6,s_tet);
2437 nod_pt->set_value(0,s_tet[0]);
2438 nod_pt->set_value(1,s_tet[1]);
2439 nod_pt->set_value(2,s_tet[2]);
2440 tet_el_pt->interpolated_x(s_tet,x_tet);
2441 nod_pt->x(0)=x_tet[0];
2442 nod_pt->x(1)=x_tet[1];
2443 nod_pt->x(2)=x_tet[2];
2446 tet_el_pt->local_coordinate_of_node(9,s_tet);
2447 nod_pt->set_value(0,s_tet[0]);
2448 nod_pt->set_value(1,s_tet[1]);
2449 nod_pt->set_value(2,s_tet[2]);
2450 tet_el_pt->interpolated_x(s_tet,x_tet);
2451 nod_pt->x(0)=x_tet[0];
2452 nod_pt->x(1)=x_tet[1];
2453 nod_pt->x(2)=x_tet[2];
2461 nod_pt->set_value(0,s_tet[0]);
2462 nod_pt->set_value(1,s_tet[1]);
2463 nod_pt->set_value(2,s_tet[2]);
2464 tet_el_pt->interpolated_x(s_tet,x_tet);
2465 nod_pt->x(0)=x_tet[0];
2466 nod_pt->x(1)=x_tet[1];
2467 nod_pt->x(2)=x_tet[2];
2470 tet_el_pt->local_coordinate_of_node(8,s_tet);
2471 nod_pt->set_value(0,s_tet[0]);
2472 nod_pt->set_value(1,s_tet[1]);
2473 nod_pt->set_value(2,s_tet[2]);
2474 tet_el_pt->interpolated_x(s_tet,x_tet);
2475 nod_pt->x(0)=x_tet[0];
2476 nod_pt->x(1)=x_tet[1];
2477 nod_pt->x(2)=x_tet[2];
2485 nod_pt->set_value(0,s_tet[0]);
2486 nod_pt->set_value(1,s_tet[1]);
2487 nod_pt->set_value(2,s_tet[2]);
2488 tet_el_pt->interpolated_x(s_tet,x_tet);
2489 nod_pt->x(0)=x_tet[0];
2490 nod_pt->x(1)=x_tet[1];
2491 nod_pt->x(2)=x_tet[2];
2499 nod_pt->set_value(0,s_tet[0]);
2500 nod_pt->set_value(1,s_tet[1]);
2501 nod_pt->set_value(2,s_tet[2]);
2502 tet_el_pt->interpolated_x(s_tet,x_tet);
2503 nod_pt->x(0)=x_tet[0];
2504 nod_pt->x(1)=x_tet[1];
2505 nod_pt->x(2)=x_tet[2];
2512 nod_pt->set_value(0,s_tet[0]);
2513 nod_pt->set_value(1,s_tet[1]);
2514 nod_pt->set_value(2,s_tet[2]);
2515 tet_el_pt->interpolated_x(s_tet,x_tet);
2516 nod_pt->x(0)=x_tet[0];
2517 nod_pt->x(1)=x_tet[1];
2518 nod_pt->x(2)=x_tet[2];
2525 FiniteElement* el_pt=
new ELEMENT;
2527 Element_pt.push_back(el_pt);
2529 TFace face0(tet_el_pt->node_pt(3),
2530 tet_el_pt->node_pt(0),
2531 tet_el_pt->node_pt(2));
2533 TFace face1(tet_el_pt->node_pt(3),
2534 tet_el_pt->node_pt(1),
2535 tet_el_pt->node_pt(2));
2537 TFace face2(tet_el_pt->node_pt(3),
2538 tet_el_pt->node_pt(1),
2539 tet_el_pt->node_pt(0));
2542 Vector<Vector<unsigned> > tet_edge_node(3);
2543 tet_edge_node[0].resize(2);
2544 tet_edge_node[0][0]=6;
2545 tet_edge_node[0][1]=0;
2546 tet_edge_node[1].resize(2);
2547 tet_edge_node[1][0]=9;
2548 tet_edge_node[1][1]=1;
2549 tet_edge_node[2].resize(2);
2550 tet_edge_node[2][0]=8;
2551 tet_edge_node[2][1]=2;
2554 unsigned central_tet_vertex=3;
2556 Node* tet_node_pt=0;
2557 Node* old_node_pt=0;
2564 tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
2565 old_node_pt=tet_node_node_pt[tet_node_pt];
2568 Node* new_node_pt=0;
2569 if (tet_node_pt->is_on_boundary())
2571 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2575 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2577 tet_node_node_pt[tet_node_pt]=new_node_pt;
2578 Node_pt.push_back(new_node_pt);
2579 Vector<double> s(3);
2580 Vector<double> s_tet(3);
2581 Vector<double> x_tet(3);
2582 el_pt->local_coordinate_of_node(j,s);
2583 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2584 tet_el_pt->interpolated_x(s_tet,x_tet);
2585 new_node_pt->x(0)=x_tet[0];
2586 new_node_pt->x(1)=x_tet[1];
2587 new_node_pt->x(2)=x_tet[2];
2592 el_pt->node_pt(j)=old_node_pt;
2602 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
2603 old_node_pt=tet_node_node_pt[tet_node_pt];
2606 Node* new_node_pt=0;
2607 if (tet_node_pt->is_on_boundary())
2609 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2613 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2615 tet_node_node_pt[tet_node_pt]=new_node_pt;
2616 Node_pt.push_back(new_node_pt);
2617 Vector<double> s(3);
2618 Vector<double> s_tet(3);
2619 Vector<double> x_tet(3);
2620 el_pt->local_coordinate_of_node(j,s);
2621 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2622 tet_el_pt->interpolated_x(s_tet,x_tet);
2623 new_node_pt->x(0)=x_tet[0];
2624 new_node_pt->x(1)=x_tet[1];
2625 new_node_pt->x(2)=x_tet[2];
2630 el_pt->node_pt(j)=old_node_pt;
2640 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
2641 old_node_pt=tet_node_node_pt[tet_node_pt];
2644 Node* new_node_pt=0;
2645 if (tet_node_pt->is_on_boundary())
2647 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2651 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2653 tet_node_node_pt[tet_node_pt]=new_node_pt;
2654 Node_pt.push_back(new_node_pt);
2655 Vector<double> s(3);
2656 Vector<double> s_tet(3);
2657 Vector<double> x_tet(3);
2658 el_pt->local_coordinate_of_node(j,s);
2659 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2660 tet_el_pt->interpolated_x(s_tet,x_tet);
2661 new_node_pt->x(0)=x_tet[0];
2662 new_node_pt->x(1)=x_tet[1];
2663 new_node_pt->x(2)=x_tet[2];
2668 el_pt->node_pt(j)=old_node_pt;
2678 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
2679 old_node_pt=tet_node_node_pt[tet_node_pt];
2682 Node* new_node_pt=0;
2683 if (tet_node_pt->is_on_boundary())
2685 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2689 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2691 tet_node_node_pt[tet_node_pt]=new_node_pt;
2692 Node_pt.push_back(new_node_pt);
2693 Vector<double> s(3);
2694 Vector<double> s_tet(3);
2695 Vector<double> x_tet(3);
2696 el_pt->local_coordinate_of_node(j,s);
2697 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2698 tet_el_pt->interpolated_x(s_tet,x_tet);
2699 new_node_pt->x(0)=x_tet[0];
2700 new_node_pt->x(1)=x_tet[1];
2701 new_node_pt->x(2)=x_tet[2];
2706 el_pt->node_pt(j)=old_node_pt;
2717 old_node_pt=tet_face_node_pt[face0];
2720 Node* new_node_pt=0;
2721 if (face0.is_boundary_face())
2723 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2727 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2729 tet_face_node_pt[face0]=new_node_pt;
2730 Node_pt.push_back(new_node_pt);
2731 Vector<double> s(3);
2732 Vector<double> s_tet(3);
2733 Vector<double> x_tet(3);
2734 el_pt->local_coordinate_of_node(j,s);
2735 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2736 tet_el_pt->interpolated_x(s_tet,x_tet);
2737 new_node_pt->x(0)=x_tet[0];
2738 new_node_pt->x(1)=x_tet[1];
2739 new_node_pt->x(2)=x_tet[2];
2744 el_pt->node_pt(j)=old_node_pt;
2753 old_node_pt=tet_face_node_pt[face1];
2756 Node* new_node_pt=0;
2757 if (face1.is_boundary_face())
2759 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2763 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2765 tet_face_node_pt[face1]=new_node_pt;
2766 Node_pt.push_back(new_node_pt);
2767 Vector<double> s(3);
2768 Vector<double> s_tet(3);
2769 Vector<double> x_tet(3);
2770 el_pt->local_coordinate_of_node(j,s);
2771 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2772 tet_el_pt->interpolated_x(s_tet,x_tet);
2773 new_node_pt->x(0)=x_tet[0];
2774 new_node_pt->x(1)=x_tet[1];
2775 new_node_pt->x(2)=x_tet[2];
2780 el_pt->node_pt(j)=old_node_pt;
2789 old_node_pt=tet_face_node_pt[face2];
2792 Node* new_node_pt=0;
2793 if (face2.is_boundary_face())
2795 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2799 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2801 tet_face_node_pt[face2]=new_node_pt;
2802 Node_pt.push_back(new_node_pt);
2803 Vector<double> s(3);
2804 Vector<double> s_tet(3);
2805 Vector<double> x_tet(3);
2806 el_pt->local_coordinate_of_node(j,s);
2807 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2808 tet_el_pt->interpolated_x(s_tet,x_tet);
2809 new_node_pt->x(0)=x_tet[0];
2810 new_node_pt->x(1)=x_tet[1];
2811 new_node_pt->x(2)=x_tet[2];
2816 el_pt->node_pt(j)=old_node_pt;
2827 el_pt->node_pt(j)=centroid_node_pt;
2837 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2838 Node_pt.push_back(new_node_pt);
2839 Vector<double> s(3);
2840 Vector<double> s_tet(3);
2841 Vector<double> x_tet(3);
2842 el_pt->local_coordinate_of_node(j,s);
2843 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2844 tet_el_pt->interpolated_x(s_tet,x_tet);
2845 new_node_pt->x(0)=x_tet[0];
2846 new_node_pt->x(1)=x_tet[1];
2847 new_node_pt->x(2)=x_tet[2];
2856 Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
2857 old_node_pt=brick_edge_node_pt[edge];
2860 Node* new_node_pt=0;
2861 if (edge.is_boundary_edge())
2863 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2867 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2869 brick_edge_node_pt[edge]=new_node_pt;
2870 Node_pt.push_back(new_node_pt);
2871 Vector<double> s(3);
2872 Vector<double> s_tet(3);
2873 Vector<double> x_tet(3);
2874 el_pt->local_coordinate_of_node(j,s);
2875 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2876 tet_el_pt->interpolated_x(s_tet,x_tet);
2877 new_node_pt->x(0)=x_tet[0];
2878 new_node_pt->x(1)=x_tet[1];
2879 new_node_pt->x(2)=x_tet[2];
2884 el_pt->node_pt(j)=old_node_pt;
2894 Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
2895 old_node_pt=brick_edge_node_pt[edge];
2898 Node* new_node_pt=0;
2899 if (edge.is_boundary_edge())
2901 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2905 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2907 brick_edge_node_pt[edge]=new_node_pt;
2908 Node_pt.push_back(new_node_pt);
2909 Vector<double> s(3);
2910 Vector<double> s_tet(3);
2911 Vector<double> x_tet(3);
2912 el_pt->local_coordinate_of_node(j,s);
2913 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2914 tet_el_pt->interpolated_x(s_tet,x_tet);
2915 new_node_pt->x(0)=x_tet[0];
2916 new_node_pt->x(1)=x_tet[1];
2917 new_node_pt->x(2)=x_tet[2];
2922 el_pt->node_pt(j)=old_node_pt;
2931 Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
2932 old_node_pt=brick_edge_node_pt[edge];
2935 Node* new_node_pt=0;
2936 if (edge.is_boundary_edge())
2938 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2942 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2944 brick_edge_node_pt[edge]=new_node_pt;
2945 Node_pt.push_back(new_node_pt);
2946 Vector<double> s(3);
2947 Vector<double> s_tet(3);
2948 Vector<double> x_tet(3);
2949 el_pt->local_coordinate_of_node(j,s);
2950 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2951 tet_el_pt->interpolated_x(s_tet,x_tet);
2952 new_node_pt->x(0)=x_tet[0];
2953 new_node_pt->x(1)=x_tet[1];
2954 new_node_pt->x(2)=x_tet[2];
2959 el_pt->node_pt(j)=old_node_pt;
2968 Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
2969 old_node_pt=brick_edge_node_pt[edge];
2972 Node* new_node_pt=0;
2973 if (edge.is_boundary_edge())
2975 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
2979 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
2981 brick_edge_node_pt[edge]=new_node_pt;
2982 Node_pt.push_back(new_node_pt);
2983 Vector<double> s(3);
2984 Vector<double> s_tet(3);
2985 Vector<double> x_tet(3);
2986 el_pt->local_coordinate_of_node(j,s);
2987 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
2988 tet_el_pt->interpolated_x(s_tet,x_tet);
2989 new_node_pt->x(0)=x_tet[0];
2990 new_node_pt->x(1)=x_tet[1];
2991 new_node_pt->x(2)=x_tet[2];
2996 el_pt->node_pt(j)=old_node_pt;
3005 Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
3006 old_node_pt=brick_edge_node_pt[edge];
3009 Node* new_node_pt=0;
3010 if (edge.is_boundary_edge())
3012 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3016 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3018 brick_edge_node_pt[edge]=new_node_pt;
3019 Node_pt.push_back(new_node_pt);
3020 Vector<double> s(3);
3021 Vector<double> s_tet(3);
3022 Vector<double> x_tet(3);
3023 el_pt->local_coordinate_of_node(j,s);
3024 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3025 tet_el_pt->interpolated_x(s_tet,x_tet);
3026 new_node_pt->x(0)=x_tet[0];
3027 new_node_pt->x(1)=x_tet[1];
3028 new_node_pt->x(2)=x_tet[2];
3033 el_pt->node_pt(j)=old_node_pt;
3043 Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
3044 old_node_pt=brick_edge_node_pt[edge];
3047 Node* new_node_pt=0;
3048 if (edge.is_boundary_edge())
3050 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3054 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3056 brick_edge_node_pt[edge]=new_node_pt;
3057 Node_pt.push_back(new_node_pt);
3058 Vector<double> s(3);
3059 Vector<double> s_tet(3);
3060 Vector<double> x_tet(3);
3061 el_pt->local_coordinate_of_node(j,s);
3062 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3063 tet_el_pt->interpolated_x(s_tet,x_tet);
3064 new_node_pt->x(0)=x_tet[0];
3065 new_node_pt->x(1)=x_tet[1];
3066 new_node_pt->x(2)=x_tet[2];
3071 el_pt->node_pt(j)=old_node_pt;
3080 Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
3081 old_node_pt=brick_edge_node_pt[edge];
3084 Node* new_node_pt=0;
3085 if (edge.is_boundary_edge())
3087 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3091 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3093 brick_edge_node_pt[edge]=new_node_pt;
3094 Node_pt.push_back(new_node_pt);
3095 Vector<double> s(3);
3096 Vector<double> s_tet(3);
3097 Vector<double> x_tet(3);
3098 el_pt->local_coordinate_of_node(j,s);
3099 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3100 tet_el_pt->interpolated_x(s_tet,x_tet);
3101 new_node_pt->x(0)=x_tet[0];
3102 new_node_pt->x(1)=x_tet[1];
3103 new_node_pt->x(2)=x_tet[2];
3108 el_pt->node_pt(j)=old_node_pt;
3118 Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
3119 old_node_pt=brick_edge_node_pt[edge];
3122 Node* new_node_pt=0;
3123 if (edge.is_boundary_edge())
3125 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3129 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3131 brick_edge_node_pt[edge]=new_node_pt;
3132 Node_pt.push_back(new_node_pt);
3133 Vector<double> s(3);
3134 Vector<double> s_tet(3);
3135 Vector<double> x_tet(3);
3136 el_pt->local_coordinate_of_node(j,s);
3137 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3138 tet_el_pt->interpolated_x(s_tet,x_tet);
3139 new_node_pt->x(0)=x_tet[0];
3140 new_node_pt->x(1)=x_tet[1];
3141 new_node_pt->x(2)=x_tet[2];
3146 el_pt->node_pt(j)=old_node_pt;
3155 Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
3156 old_node_pt=brick_edge_node_pt[edge];
3159 Node* new_node_pt=0;
3160 if (edge.is_boundary_edge())
3162 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3166 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3168 brick_edge_node_pt[edge]=new_node_pt;
3169 Node_pt.push_back(new_node_pt);
3170 Vector<double> s(3);
3171 Vector<double> s_tet(3);
3172 Vector<double> x_tet(3);
3173 el_pt->local_coordinate_of_node(j,s);
3174 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3175 tet_el_pt->interpolated_x(s_tet,x_tet);
3176 new_node_pt->x(0)=x_tet[0];
3177 new_node_pt->x(1)=x_tet[1];
3178 new_node_pt->x(2)=x_tet[2];
3183 el_pt->node_pt(j)=old_node_pt;
3193 Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
3194 old_node_pt=brick_edge_node_pt[edge];
3197 Node* new_node_pt=0;
3198 if (edge.is_boundary_edge())
3200 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3204 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3206 brick_edge_node_pt[edge]=new_node_pt;
3207 Node_pt.push_back(new_node_pt);
3208 Vector<double> s(3);
3209 Vector<double> s_tet(3);
3210 Vector<double> x_tet(3);
3211 el_pt->local_coordinate_of_node(j,s);
3212 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3213 tet_el_pt->interpolated_x(s_tet,x_tet);
3214 new_node_pt->x(0)=x_tet[0];
3215 new_node_pt->x(1)=x_tet[1];
3216 new_node_pt->x(2)=x_tet[2];
3221 el_pt->node_pt(j)=old_node_pt;
3231 Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
3232 old_node_pt=brick_edge_node_pt[edge];
3235 Node* new_node_pt=0;
3236 if (edge.is_boundary_edge())
3238 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3242 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3244 brick_edge_node_pt[edge]=new_node_pt;
3245 Node_pt.push_back(new_node_pt);
3246 Vector<double> s(3);
3247 Vector<double> s_tet(3);
3248 Vector<double> x_tet(3);
3249 el_pt->local_coordinate_of_node(j,s);
3250 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3251 tet_el_pt->interpolated_x(s_tet,x_tet);
3252 new_node_pt->x(0)=x_tet[0];
3253 new_node_pt->x(1)=x_tet[1];
3254 new_node_pt->x(2)=x_tet[2];
3259 el_pt->node_pt(j)=old_node_pt;
3269 Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
3270 old_node_pt=brick_edge_node_pt[edge];
3273 Node* new_node_pt=0;
3274 if (edge.is_boundary_edge())
3276 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3280 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3282 brick_edge_node_pt[edge]=new_node_pt;
3283 Node_pt.push_back(new_node_pt);
3284 Vector<double> s(3);
3285 Vector<double> s_tet(3);
3286 Vector<double> x_tet(3);
3287 el_pt->local_coordinate_of_node(j,s);
3288 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3289 tet_el_pt->interpolated_x(s_tet,x_tet);
3290 new_node_pt->x(0)=x_tet[0];
3291 new_node_pt->x(1)=x_tet[1];
3292 new_node_pt->x(2)=x_tet[2];
3297 el_pt->node_pt(j)=old_node_pt;
3308 TFace face(tet_el_pt->node_pt( central_tet_vertex),
3309 tet_el_pt->node_pt(tet_edge_node[0][0]),
3310 tet_el_pt->node_pt(tet_edge_node[2][0]));
3312 old_node_pt=tet_face_node_pt[face];
3315 Node* new_node_pt=0;
3316 if (face.is_boundary_face())
3318 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3322 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3324 tet_face_node_pt[face]=new_node_pt;
3325 Node_pt.push_back(new_node_pt);
3326 Vector<double> s(3);
3327 Vector<double> s_tet(3);
3328 Vector<double> x_tet(3);
3329 el_pt->local_coordinate_of_node(j,s);
3330 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3331 tet_el_pt->interpolated_x(s_tet,x_tet);
3332 new_node_pt->x(0)=x_tet[0];
3333 new_node_pt->x(1)=x_tet[1];
3334 new_node_pt->x(2)=x_tet[2];
3339 el_pt->node_pt(j)=old_node_pt;
3351 TFace face(tet_el_pt->node_pt(central_tet_vertex),
3352 tet_el_pt->node_pt(tet_edge_node[1][0]),
3353 tet_el_pt->node_pt(tet_edge_node[2][0]));
3354 old_node_pt=tet_face_node_pt[face];
3357 Node* new_node_pt=0;
3358 if (face.is_boundary_face())
3360 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3364 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3366 tet_face_node_pt[face]=new_node_pt;
3367 Node_pt.push_back(new_node_pt);
3368 Vector<double> s(3);
3369 Vector<double> s_tet(3);
3370 Vector<double> x_tet(3);
3371 el_pt->local_coordinate_of_node(j,s);
3372 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3373 tet_el_pt->interpolated_x(s_tet,x_tet);
3374 new_node_pt->x(0)=x_tet[0];
3375 new_node_pt->x(1)=x_tet[1];
3376 new_node_pt->x(2)=x_tet[2];
3381 el_pt->node_pt(j)=old_node_pt;
3393 TFace face(tet_el_pt->node_pt(central_tet_vertex),
3394 tet_el_pt->node_pt(tet_edge_node[0][0]),
3395 tet_el_pt->node_pt(tet_edge_node[1][0]));
3397 old_node_pt=tet_face_node_pt[face];
3400 Node* new_node_pt=0;
3401 if (face.is_boundary_face())
3403 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3407 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3409 tet_face_node_pt[face]=new_node_pt;
3410 Node_pt.push_back(new_node_pt);
3411 Vector<double> s(3);
3412 Vector<double> s_tet(3);
3413 Vector<double> x_tet(3);
3414 el_pt->local_coordinate_of_node(j,s);
3415 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3416 tet_el_pt->interpolated_x(s_tet,x_tet);
3417 new_node_pt->x(0)=x_tet[0];
3418 new_node_pt->x(1)=x_tet[1];
3419 new_node_pt->x(2)=x_tet[2];
3424 el_pt->node_pt(j)=old_node_pt;
3432 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3433 Node_pt.push_back(new_node_pt);
3434 Vector<double> s(3);
3435 Vector<double> s_tet(3);
3436 Vector<double> x_tet(3);
3437 el_pt->local_coordinate_of_node(j,s);
3438 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
3439 top_mid_face_node2_pt=new_node_pt;
3440 tet_el_pt->interpolated_x(s_tet,x_tet);
3441 new_node_pt->x(0)=x_tet[0];
3442 new_node_pt->x(1)=x_tet[1];
3443 new_node_pt->x(2)=x_tet[2];
3452 el_pt->node_pt(j)=back_mid_face_node0_pt;
3461 el_pt->node_pt(j)=right_mid_face_node1_pt;
3471 for (
unsigned j=0;j<8;j++)
3473 Node* nod_pt=dummy_q_el_pt[3]->node_pt(j);
3474 Vector<double> s_tet(3);
3475 Vector<double> x_tet(3);
3479 tet_el_pt->local_coordinate_of_node(2,s_tet);
3480 nod_pt->set_value(0,s_tet[0]);
3481 nod_pt->set_value(1,s_tet[1]);
3482 nod_pt->set_value(2,s_tet[2]);
3483 tet_el_pt->interpolated_x(s_tet,x_tet);
3484 nod_pt->x(0)=x_tet[0];
3485 nod_pt->x(1)=x_tet[1];
3486 nod_pt->x(2)=x_tet[2];
3489 tet_el_pt->local_coordinate_of_node(7,s_tet);
3490 nod_pt->set_value(0,s_tet[0]);
3491 nod_pt->set_value(1,s_tet[1]);
3492 nod_pt->set_value(2,s_tet[2]);
3493 tet_el_pt->interpolated_x(s_tet,x_tet);
3494 nod_pt->x(0)=x_tet[0];
3495 nod_pt->x(1)=x_tet[1];
3496 nod_pt->x(2)=x_tet[2];
3499 tet_el_pt->local_coordinate_of_node(5,s_tet);
3500 nod_pt->set_value(0,s_tet[0]);
3501 nod_pt->set_value(1,s_tet[1]);
3502 nod_pt->set_value(2,s_tet[2]);
3503 tet_el_pt->interpolated_x(s_tet,x_tet);
3504 nod_pt->x(0)=x_tet[0];
3505 nod_pt->x(1)=x_tet[1];
3506 nod_pt->x(2)=x_tet[2];
3514 nod_pt->set_value(0,s_tet[0]);
3515 nod_pt->set_value(1,s_tet[1]);
3516 nod_pt->set_value(2,s_tet[2]);
3517 tet_el_pt->interpolated_x(s_tet,x_tet);
3518 nod_pt->x(0)=x_tet[0];
3519 nod_pt->x(1)=x_tet[1];
3520 nod_pt->x(2)=x_tet[2];
3523 tet_el_pt->local_coordinate_of_node(8,s_tet);
3524 nod_pt->set_value(0,s_tet[0]);
3525 nod_pt->set_value(1,s_tet[1]);
3526 nod_pt->set_value(2,s_tet[2]);
3527 tet_el_pt->interpolated_x(s_tet,x_tet);
3528 nod_pt->x(0)=x_tet[0];
3529 nod_pt->x(1)=x_tet[1];
3530 nod_pt->x(2)=x_tet[2];
3538 nod_pt->set_value(0,s_tet[0]);
3539 nod_pt->set_value(1,s_tet[1]);
3540 nod_pt->set_value(2,s_tet[2]);
3541 tet_el_pt->interpolated_x(s_tet,x_tet);
3542 nod_pt->x(0)=x_tet[0];
3543 nod_pt->x(1)=x_tet[1];
3544 nod_pt->x(2)=x_tet[2];
3552 nod_pt->set_value(0,s_tet[0]);
3553 nod_pt->set_value(1,s_tet[1]);
3554 nod_pt->set_value(2,s_tet[2]);
3555 tet_el_pt->interpolated_x(s_tet,x_tet);
3556 nod_pt->x(0)=x_tet[0];
3557 nod_pt->x(1)=x_tet[1];
3558 nod_pt->x(2)=x_tet[2];
3565 nod_pt->set_value(0,s_tet[0]);
3566 nod_pt->set_value(1,s_tet[1]);
3567 nod_pt->set_value(2,s_tet[2]);
3568 tet_el_pt->interpolated_x(s_tet,x_tet);
3569 nod_pt->x(0)=x_tet[0];
3570 nod_pt->x(1)=x_tet[1];
3571 nod_pt->x(2)=x_tet[2];
3579 FiniteElement* el_pt=
new ELEMENT;
3581 Element_pt.push_back(el_pt);
3583 TFace face0(tet_el_pt->node_pt(1),
3584 tet_el_pt->node_pt(2),
3585 tet_el_pt->node_pt(3));
3587 TFace face1(tet_el_pt->node_pt(0),
3588 tet_el_pt->node_pt(2),
3589 tet_el_pt->node_pt(3));
3591 TFace face2(tet_el_pt->node_pt(0),
3592 tet_el_pt->node_pt(1),
3593 tet_el_pt->node_pt(2));
3596 Vector<Vector<unsigned> > tet_edge_node(3);
3597 tet_edge_node[0].resize(2);
3598 tet_edge_node[0][0]=7;
3599 tet_edge_node[0][1]=1;
3600 tet_edge_node[1].resize(2);
3601 tet_edge_node[1][0]=5;
3602 tet_edge_node[1][1]=0;
3603 tet_edge_node[2].resize(2);
3604 tet_edge_node[2][0]=8;
3605 tet_edge_node[2][1]=3;
3608 unsigned central_tet_vertex=2;
3610 Node* tet_node_pt=0;
3611 Node* old_node_pt=0;
3618 tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
3619 old_node_pt=tet_node_node_pt[tet_node_pt];
3622 Node* new_node_pt=0;
3623 if (tet_node_pt->is_on_boundary())
3625 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3629 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3631 tet_node_node_pt[tet_node_pt]=new_node_pt;
3632 Node_pt.push_back(new_node_pt);
3633 Vector<double> s(3);
3634 Vector<double> s_tet(3);
3635 Vector<double> x_tet(3);
3636 el_pt->local_coordinate_of_node(j,s);
3637 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3638 tet_el_pt->interpolated_x(s_tet,x_tet);
3639 new_node_pt->x(0)=x_tet[0];
3640 new_node_pt->x(1)=x_tet[1];
3641 new_node_pt->x(2)=x_tet[2];
3646 el_pt->node_pt(j)=old_node_pt;
3656 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
3657 old_node_pt=tet_node_node_pt[tet_node_pt];
3660 Node* new_node_pt=0;
3661 if (tet_node_pt->is_on_boundary())
3663 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3667 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3669 tet_node_node_pt[tet_node_pt]=new_node_pt;
3670 Node_pt.push_back(new_node_pt);
3671 Vector<double> s(3);
3672 Vector<double> s_tet(3);
3673 Vector<double> x_tet(3);
3674 el_pt->local_coordinate_of_node(j,s);
3675 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3676 tet_el_pt->interpolated_x(s_tet,x_tet);
3677 new_node_pt->x(0)=x_tet[0];
3678 new_node_pt->x(1)=x_tet[1];
3679 new_node_pt->x(2)=x_tet[2];
3684 el_pt->node_pt(j)=old_node_pt;
3694 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
3695 old_node_pt=tet_node_node_pt[tet_node_pt];
3698 Node* new_node_pt=0;
3699 if (tet_node_pt->is_on_boundary())
3701 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3705 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3707 tet_node_node_pt[tet_node_pt]=new_node_pt;
3708 Node_pt.push_back(new_node_pt);
3709 Vector<double> s(3);
3710 Vector<double> s_tet(3);
3711 Vector<double> x_tet(3);
3712 el_pt->local_coordinate_of_node(j,s);
3713 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3714 tet_el_pt->interpolated_x(s_tet,x_tet);
3715 new_node_pt->x(0)=x_tet[0];
3716 new_node_pt->x(1)=x_tet[1];
3717 new_node_pt->x(2)=x_tet[2];
3722 el_pt->node_pt(j)=old_node_pt;
3732 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
3733 old_node_pt=tet_node_node_pt[tet_node_pt];
3736 Node* new_node_pt=0;
3737 if (tet_node_pt->is_on_boundary())
3739 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3743 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3745 tet_node_node_pt[tet_node_pt]=new_node_pt;
3746 Node_pt.push_back(new_node_pt);
3747 Vector<double> s(3);
3748 Vector<double> s_tet(3);
3749 Vector<double> x_tet(3);
3750 el_pt->local_coordinate_of_node(j,s);
3751 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3752 tet_el_pt->interpolated_x(s_tet,x_tet);
3753 new_node_pt->x(0)=x_tet[0];
3754 new_node_pt->x(1)=x_tet[1];
3755 new_node_pt->x(2)=x_tet[2];
3760 el_pt->node_pt(j)=old_node_pt;
3771 old_node_pt=tet_face_node_pt[face0];
3774 Node* new_node_pt=0;
3775 if (face0.is_boundary_face())
3777 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3781 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3783 tet_face_node_pt[face0]=new_node_pt;
3784 Node_pt.push_back(new_node_pt);
3785 Vector<double> s(3);
3786 Vector<double> s_tet(3);
3787 Vector<double> x_tet(3);
3788 el_pt->local_coordinate_of_node(j,s);
3789 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3790 tet_el_pt->interpolated_x(s_tet,x_tet);
3791 new_node_pt->x(0)=x_tet[0];
3792 new_node_pt->x(1)=x_tet[1];
3793 new_node_pt->x(2)=x_tet[2];
3798 el_pt->node_pt(j)=old_node_pt;
3807 old_node_pt=tet_face_node_pt[face1];
3810 Node* new_node_pt=0;
3811 if (face1.is_boundary_face())
3813 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3817 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3819 tet_face_node_pt[face1]=new_node_pt;
3820 Node_pt.push_back(new_node_pt);
3821 Vector<double> s(3);
3822 Vector<double> s_tet(3);
3823 Vector<double> x_tet(3);
3824 el_pt->local_coordinate_of_node(j,s);
3825 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3826 tet_el_pt->interpolated_x(s_tet,x_tet);
3827 new_node_pt->x(0)=x_tet[0];
3828 new_node_pt->x(1)=x_tet[1];
3829 new_node_pt->x(2)=x_tet[2];
3834 el_pt->node_pt(j)=old_node_pt;
3843 old_node_pt=tet_face_node_pt[face2];
3846 Node* new_node_pt=0;
3847 if (face2.is_boundary_face())
3849 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3853 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3855 tet_face_node_pt[face2]=new_node_pt;
3856 Node_pt.push_back(new_node_pt);
3857 Vector<double> s(3);
3858 Vector<double> s_tet(3);
3859 Vector<double> x_tet(3);
3860 el_pt->local_coordinate_of_node(j,s);
3861 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3862 tet_el_pt->interpolated_x(s_tet,x_tet);
3863 new_node_pt->x(0)=x_tet[0];
3864 new_node_pt->x(1)=x_tet[1];
3865 new_node_pt->x(2)=x_tet[2];
3870 el_pt->node_pt(j)=old_node_pt;
3881 el_pt->node_pt(j)=centroid_node_pt;
3891 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3892 Node_pt.push_back(new_node_pt);
3893 Vector<double> s(3);
3894 Vector<double> s_tet(3);
3895 Vector<double> x_tet(3);
3896 el_pt->local_coordinate_of_node(j,s);
3897 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3898 tet_el_pt->interpolated_x(s_tet,x_tet);
3899 new_node_pt->x(0)=x_tet[0];
3900 new_node_pt->x(1)=x_tet[1];
3901 new_node_pt->x(2)=x_tet[2];
3910 Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
3911 old_node_pt=brick_edge_node_pt[edge];
3914 Node* new_node_pt=0;
3915 if (edge.is_boundary_edge())
3917 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3921 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3923 brick_edge_node_pt[edge]=new_node_pt;
3924 Node_pt.push_back(new_node_pt);
3925 Vector<double> s(3);
3926 Vector<double> s_tet(3);
3927 Vector<double> x_tet(3);
3928 el_pt->local_coordinate_of_node(j,s);
3929 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3930 tet_el_pt->interpolated_x(s_tet,x_tet);
3931 new_node_pt->x(0)=x_tet[0];
3932 new_node_pt->x(1)=x_tet[1];
3933 new_node_pt->x(2)=x_tet[2];
3938 el_pt->node_pt(j)=old_node_pt;
3948 Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
3949 old_node_pt=brick_edge_node_pt[edge];
3952 Node* new_node_pt=0;
3953 if (edge.is_boundary_edge())
3955 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3959 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3961 brick_edge_node_pt[edge]=new_node_pt;
3962 Node_pt.push_back(new_node_pt);
3963 Vector<double> s(3);
3964 Vector<double> s_tet(3);
3965 Vector<double> x_tet(3);
3966 el_pt->local_coordinate_of_node(j,s);
3967 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
3968 tet_el_pt->interpolated_x(s_tet,x_tet);
3969 new_node_pt->x(0)=x_tet[0];
3970 new_node_pt->x(1)=x_tet[1];
3971 new_node_pt->x(2)=x_tet[2];
3976 el_pt->node_pt(j)=old_node_pt;
3985 Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
3986 old_node_pt=brick_edge_node_pt[edge];
3989 Node* new_node_pt=0;
3990 if (edge.is_boundary_edge())
3992 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
3996 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
3998 brick_edge_node_pt[edge]=new_node_pt;
3999 Node_pt.push_back(new_node_pt);
4000 Vector<double> s(3);
4001 Vector<double> s_tet(3);
4002 Vector<double> x_tet(3);
4003 el_pt->local_coordinate_of_node(j,s);
4004 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4005 tet_el_pt->interpolated_x(s_tet,x_tet);
4006 new_node_pt->x(0)=x_tet[0];
4007 new_node_pt->x(1)=x_tet[1];
4008 new_node_pt->x(2)=x_tet[2];
4013 el_pt->node_pt(j)=old_node_pt;
4022 Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
4023 old_node_pt=brick_edge_node_pt[edge];
4026 Node* new_node_pt=0;
4027 if (edge.is_boundary_edge())
4029 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4033 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4035 brick_edge_node_pt[edge]=new_node_pt;
4036 Node_pt.push_back(new_node_pt);
4037 Vector<double> s(3);
4038 Vector<double> s_tet(3);
4039 Vector<double> x_tet(3);
4040 el_pt->local_coordinate_of_node(j,s);
4041 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4042 tet_el_pt->interpolated_x(s_tet,x_tet);
4043 new_node_pt->x(0)=x_tet[0];
4044 new_node_pt->x(1)=x_tet[1];
4045 new_node_pt->x(2)=x_tet[2];
4050 el_pt->node_pt(j)=old_node_pt;
4059 Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
4060 old_node_pt=brick_edge_node_pt[edge];
4063 Node* new_node_pt=0;
4064 if (edge.is_boundary_edge())
4066 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4070 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4072 brick_edge_node_pt[edge]=new_node_pt;
4073 Node_pt.push_back(new_node_pt);
4074 Vector<double> s(3);
4075 Vector<double> s_tet(3);
4076 Vector<double> x_tet(3);
4077 el_pt->local_coordinate_of_node(j,s);
4078 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4079 tet_el_pt->interpolated_x(s_tet,x_tet);
4080 new_node_pt->x(0)=x_tet[0];
4081 new_node_pt->x(1)=x_tet[1];
4082 new_node_pt->x(2)=x_tet[2];
4087 el_pt->node_pt(j)=old_node_pt;
4097 Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
4098 old_node_pt=brick_edge_node_pt[edge];
4101 Node* new_node_pt=0;
4102 if (edge.is_boundary_edge())
4104 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4108 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4110 brick_edge_node_pt[edge]=new_node_pt;
4111 Node_pt.push_back(new_node_pt);
4112 Vector<double> s(3);
4113 Vector<double> s_tet(3);
4114 Vector<double> x_tet(3);
4115 el_pt->local_coordinate_of_node(j,s);
4116 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4117 tet_el_pt->interpolated_x(s_tet,x_tet);
4118 new_node_pt->x(0)=x_tet[0];
4119 new_node_pt->x(1)=x_tet[1];
4120 new_node_pt->x(2)=x_tet[2];
4125 el_pt->node_pt(j)=old_node_pt;
4134 Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
4135 old_node_pt=brick_edge_node_pt[edge];
4138 Node* new_node_pt=0;
4139 if (edge.is_boundary_edge())
4141 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4145 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4147 brick_edge_node_pt[edge]=new_node_pt;
4148 Node_pt.push_back(new_node_pt);
4149 Vector<double> s(3);
4150 Vector<double> s_tet(3);
4151 Vector<double> x_tet(3);
4152 el_pt->local_coordinate_of_node(j,s);
4153 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4154 tet_el_pt->interpolated_x(s_tet,x_tet);
4155 new_node_pt->x(0)=x_tet[0];
4156 new_node_pt->x(1)=x_tet[1];
4157 new_node_pt->x(2)=x_tet[2];
4162 el_pt->node_pt(j)=old_node_pt;
4172 Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
4173 old_node_pt=brick_edge_node_pt[edge];
4176 Node* new_node_pt=0;
4177 if (edge.is_boundary_edge())
4179 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4183 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4185 brick_edge_node_pt[edge]=new_node_pt;
4186 Node_pt.push_back(new_node_pt);
4187 Vector<double> s(3);
4188 Vector<double> s_tet(3);
4189 Vector<double> x_tet(3);
4190 el_pt->local_coordinate_of_node(j,s);
4191 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4192 tet_el_pt->interpolated_x(s_tet,x_tet);
4193 new_node_pt->x(0)=x_tet[0];
4194 new_node_pt->x(1)=x_tet[1];
4195 new_node_pt->x(2)=x_tet[2];
4200 el_pt->node_pt(j)=old_node_pt;
4209 Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
4210 old_node_pt=brick_edge_node_pt[edge];
4213 Node* new_node_pt=0;
4214 if (edge.is_boundary_edge())
4216 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4220 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4222 brick_edge_node_pt[edge]=new_node_pt;
4223 Node_pt.push_back(new_node_pt);
4224 Vector<double> s(3);
4225 Vector<double> s_tet(3);
4226 Vector<double> x_tet(3);
4227 el_pt->local_coordinate_of_node(j,s);
4228 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4229 tet_el_pt->interpolated_x(s_tet,x_tet);
4230 new_node_pt->x(0)=x_tet[0];
4231 new_node_pt->x(1)=x_tet[1];
4232 new_node_pt->x(2)=x_tet[2];
4237 el_pt->node_pt(j)=old_node_pt;
4247 Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
4248 old_node_pt=brick_edge_node_pt[edge];
4251 Node* new_node_pt=0;
4252 if (edge.is_boundary_edge())
4254 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4258 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4260 brick_edge_node_pt[edge]=new_node_pt;
4261 Node_pt.push_back(new_node_pt);
4262 Vector<double> s(3);
4263 Vector<double> s_tet(3);
4264 Vector<double> x_tet(3);
4265 el_pt->local_coordinate_of_node(j,s);
4266 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4267 tet_el_pt->interpolated_x(s_tet,x_tet);
4268 new_node_pt->x(0)=x_tet[0];
4269 new_node_pt->x(1)=x_tet[1];
4270 new_node_pt->x(2)=x_tet[2];
4275 el_pt->node_pt(j)=old_node_pt;
4285 Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
4286 old_node_pt=brick_edge_node_pt[edge];
4289 Node* new_node_pt=0;
4290 if (edge.is_boundary_edge())
4292 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4296 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4298 brick_edge_node_pt[edge]=new_node_pt;
4299 Node_pt.push_back(new_node_pt);
4300 Vector<double> s(3);
4301 Vector<double> s_tet(3);
4302 Vector<double> x_tet(3);
4303 el_pt->local_coordinate_of_node(j,s);
4304 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4305 tet_el_pt->interpolated_x(s_tet,x_tet);
4306 new_node_pt->x(0)=x_tet[0];
4307 new_node_pt->x(1)=x_tet[1];
4308 new_node_pt->x(2)=x_tet[2];
4313 el_pt->node_pt(j)=old_node_pt;
4323 Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
4324 old_node_pt=brick_edge_node_pt[edge];
4327 Node* new_node_pt=0;
4328 if (edge.is_boundary_edge())
4330 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4334 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4336 brick_edge_node_pt[edge]=new_node_pt;
4337 Node_pt.push_back(new_node_pt);
4338 Vector<double> s(3);
4339 Vector<double> s_tet(3);
4340 Vector<double> x_tet(3);
4341 el_pt->local_coordinate_of_node(j,s);
4342 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4343 tet_el_pt->interpolated_x(s_tet,x_tet);
4344 new_node_pt->x(0)=x_tet[0];
4345 new_node_pt->x(1)=x_tet[1];
4346 new_node_pt->x(2)=x_tet[2];
4351 el_pt->node_pt(j)=old_node_pt;
4362 TFace face(tet_el_pt->node_pt( central_tet_vertex),
4363 tet_el_pt->node_pt(tet_edge_node[0][0]),
4364 tet_el_pt->node_pt(tet_edge_node[2][0]));
4366 old_node_pt=tet_face_node_pt[face];
4369 Node* new_node_pt=0;
4370 if (face.is_boundary_face())
4372 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4376 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4378 tet_face_node_pt[face]=new_node_pt;
4379 Node_pt.push_back(new_node_pt);
4380 Vector<double> s(3);
4381 Vector<double> s_tet(3);
4382 Vector<double> x_tet(3);
4383 el_pt->local_coordinate_of_node(j,s);
4384 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4385 tet_el_pt->interpolated_x(s_tet,x_tet);
4386 new_node_pt->x(0)=x_tet[0];
4387 new_node_pt->x(1)=x_tet[1];
4388 new_node_pt->x(2)=x_tet[2];
4393 el_pt->node_pt(j)=old_node_pt;
4405 TFace face(tet_el_pt->node_pt(central_tet_vertex),
4406 tet_el_pt->node_pt(tet_edge_node[1][0]),
4407 tet_el_pt->node_pt(tet_edge_node[2][0]));
4409 old_node_pt=tet_face_node_pt[face];
4412 Node* new_node_pt=0;
4413 if (face.is_boundary_face())
4415 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4419 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4421 tet_face_node_pt[face]=new_node_pt;
4422 Node_pt.push_back(new_node_pt);
4423 Vector<double> s(3);
4424 Vector<double> s_tet(3);
4425 Vector<double> x_tet(3);
4426 el_pt->local_coordinate_of_node(j,s);
4427 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4428 tet_el_pt->interpolated_x(s_tet,x_tet);
4429 new_node_pt->x(0)=x_tet[0];
4430 new_node_pt->x(1)=x_tet[1];
4431 new_node_pt->x(2)=x_tet[2];
4436 el_pt->node_pt(j)=old_node_pt;
4448 TFace face(tet_el_pt->node_pt(central_tet_vertex),
4449 tet_el_pt->node_pt(tet_edge_node[0][0]),
4450 tet_el_pt->node_pt(tet_edge_node[1][0]));
4452 old_node_pt=tet_face_node_pt[face];
4455 Node* new_node_pt=0;
4456 if (face.is_boundary_face())
4458 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
4462 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
4464 tet_face_node_pt[face]=new_node_pt;
4465 Node_pt.push_back(new_node_pt);
4466 Vector<double> s(3);
4467 Vector<double> s_tet(3);
4468 Vector<double> x_tet(3);
4469 el_pt->local_coordinate_of_node(j,s);
4470 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
4471 tet_el_pt->interpolated_x(s_tet,x_tet);
4472 new_node_pt->x(0)=x_tet[0];
4473 new_node_pt->x(1)=x_tet[1];
4474 new_node_pt->x(2)=x_tet[2];
4479 el_pt->node_pt(j)=old_node_pt;
4489 el_pt->node_pt(j)=top_mid_face_node2_pt;
4498 el_pt->node_pt(j)=top_mid_face_node1_pt;
4507 el_pt->node_pt(j)=top_mid_face_node0_pt;
4519 for (
int face_index=0;face_index<4;face_index++)
4529 face_pt=
new TFace(tet_el_pt->node_pt(1),
4530 tet_el_pt->node_pt(2),
4531 tet_el_pt->node_pt(3));
4536 face_pt=
new TFace(tet_el_pt->node_pt(0),
4537 tet_el_pt->node_pt(2),
4538 tet_el_pt->node_pt(3));
4545 face_pt=
new TFace(tet_el_pt->node_pt(0),
4546 tet_el_pt->node_pt(1),
4547 tet_el_pt->node_pt(3));
4553 face_pt=
new TFace(tet_el_pt->node_pt(0),
4554 tet_el_pt->node_pt(1),
4555 tet_el_pt->node_pt(2));
4561 if (face_pt->is_boundary_face())
4563 std::set<unsigned> bnd;
4564 std::set<unsigned>* bnd_pt=&bnd;
4565 face_pt->get_boundaries_pt(bnd_pt);
4568 if ((*bnd_pt).size()>1)
4570 std::ostringstream error_stream;
4571 error_stream <<
"TFace should only be on one boundary.\n";
4572 throw OomphLibError(
4574 OOMPH_CURRENT_FUNCTION,
4575 OOMPH_EXCEPTION_LOCATION);
4579 if ((*bnd_pt).size()==1)
4582 FaceElement* face_el_pt=0;
4583 if (tet_mesh_is_solid_mesh)
4586 if (
dynamic_cast<SolidTElement<3,3>*
>(tet_el_pt)==0)
4588 std::ostringstream error_stream;
4590 <<
"Tet-element cannot be cast to SolidTElement<3,3>.\n" 4591 <<
"BrickFromTetMesh can only be built from\n" 4592 <<
"mesh containing quadratic tets.\n" 4594 throw OomphLibError(error_stream.str(),
4595 OOMPH_CURRENT_FUNCTION,
4596 OOMPH_EXCEPTION_LOCATION);
4601 new DummyFaceElement<SolidTElement<3,3> >(tet_el_pt,face_index);
4606 if (
dynamic_cast<TElement<3,3>*
>(tet_el_pt)==0)
4608 std::ostringstream error_stream;
4610 <<
"Tet-element cannot be cast to TElement<3,3>.\n" 4611 <<
"BrickFromTetMesh can only be built from\n" 4612 <<
"mesh containing quadratic tets.\n" 4614 throw OomphLibError(error_stream.str(),
4615 OOMPH_CURRENT_FUNCTION,
4616 OOMPH_EXCEPTION_LOCATION);
4622 new DummyFaceElement<TElement<3,3> >(tet_el_pt,face_index);
4628 unsigned b=(*(*bnd_pt).begin());
4629 Boundary_coordinate_exists[b]=
true;
4630 face_el_pt->set_boundary_number_in_bulk_mesh(b);
4634 Vector<Node*> brick_face_node_pt(19);
4639 brick_face_node_pt[0]=brick_el1_pt->node_pt(0);
4640 brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
4641 brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
4643 brick_face_node_pt[3]=brick_el1_pt->node_pt(18);
4644 brick_face_node_pt[4]=brick_el2_pt->node_pt(18);
4645 brick_face_node_pt[5]=brick_el1_pt->node_pt(2);
4647 brick_face_node_pt[6]=brick_el1_pt->node_pt(9);
4648 brick_face_node_pt[7]=brick_el3_pt->node_pt(1);
4649 brick_face_node_pt[8]=brick_el3_pt->node_pt(9);
4651 brick_face_node_pt[9]=brick_el2_pt->node_pt(9);
4652 brick_face_node_pt[10]=brick_el2_pt->node_pt(3);
4653 brick_face_node_pt[11]=brick_el1_pt->node_pt(1);
4655 brick_face_node_pt[12]=brick_el1_pt->node_pt(20);
4657 brick_face_node_pt[13]=brick_el2_pt->node_pt(12);
4658 brick_face_node_pt[14]=brick_el1_pt->node_pt(19);
4660 brick_face_node_pt[15]=brick_el1_pt->node_pt(10);
4661 brick_face_node_pt[16]=brick_el2_pt->node_pt(21);
4663 brick_face_node_pt[17]=brick_el3_pt->node_pt(10);
4664 brick_face_node_pt[18]=brick_el1_pt->node_pt(11);
4668 brick_face_node_pt[0]=brick_el0_pt->node_pt(0);
4669 brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
4670 brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
4672 brick_face_node_pt[3]=brick_el0_pt->node_pt(18);
4673 brick_face_node_pt[4]=brick_el2_pt->node_pt(18);
4674 brick_face_node_pt[5]=brick_el0_pt->node_pt(6);
4676 brick_face_node_pt[6]=brick_el0_pt->node_pt(9);
4677 brick_face_node_pt[7]=brick_el3_pt->node_pt(3);
4678 brick_face_node_pt[8]=brick_el3_pt->node_pt(9);
4680 brick_face_node_pt[9]=brick_el2_pt->node_pt(9);
4681 brick_face_node_pt[10]=brick_el2_pt->node_pt(1);
4682 brick_face_node_pt[11]=brick_el0_pt->node_pt(3);
4684 brick_face_node_pt[12]=brick_el0_pt->node_pt(24);
4686 brick_face_node_pt[13]=brick_el2_pt->node_pt(10);
4687 brick_face_node_pt[14]=brick_el0_pt->node_pt(21);
4689 brick_face_node_pt[15]=brick_el0_pt->node_pt(12);
4690 brick_face_node_pt[16]=brick_el3_pt->node_pt(21);
4692 brick_face_node_pt[17]=brick_el3_pt->node_pt(12);
4693 brick_face_node_pt[18]=brick_el0_pt->node_pt(15);
4697 brick_face_node_pt[0]=brick_el0_pt->node_pt(0);
4698 brick_face_node_pt[1]=brick_el1_pt->node_pt(0);
4699 brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
4701 brick_face_node_pt[3]=brick_el0_pt->node_pt(2);
4702 brick_face_node_pt[4]=brick_el1_pt->node_pt(2);
4703 brick_face_node_pt[5]=brick_el0_pt->node_pt(6);
4705 brick_face_node_pt[6]=brick_el0_pt->node_pt(1);
4706 brick_face_node_pt[7]=brick_el1_pt->node_pt(3);
4707 brick_face_node_pt[8]=brick_el1_pt->node_pt(1);
4709 brick_face_node_pt[9]=brick_el2_pt->node_pt(3);
4710 brick_face_node_pt[10]=brick_el2_pt->node_pt(1);
4711 brick_face_node_pt[11]=brick_el0_pt->node_pt(3);
4713 brick_face_node_pt[12]=brick_el0_pt->node_pt(8);
4715 brick_face_node_pt[13]=brick_el2_pt->node_pt(4);
4716 brick_face_node_pt[14]=brick_el0_pt->node_pt(5);
4718 brick_face_node_pt[15]=brick_el0_pt->node_pt(4);
4719 brick_face_node_pt[16]=brick_el1_pt->node_pt(5);
4721 brick_face_node_pt[17]=brick_el1_pt->node_pt(4);
4722 brick_face_node_pt[18]=brick_el0_pt->node_pt(7);
4726 brick_face_node_pt[0]=brick_el1_pt->node_pt(0);
4727 brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
4728 brick_face_node_pt[2]=brick_el0_pt->node_pt(0);
4730 brick_face_node_pt[3]=brick_el1_pt->node_pt(18);
4731 brick_face_node_pt[4]=brick_el3_pt->node_pt(6);
4732 brick_face_node_pt[5]=brick_el1_pt->node_pt(6);
4734 brick_face_node_pt[6]=brick_el1_pt->node_pt(9);
4735 brick_face_node_pt[7]=brick_el3_pt->node_pt(1);
4736 brick_face_node_pt[8]=brick_el3_pt->node_pt(3);
4738 brick_face_node_pt[9]=brick_el0_pt->node_pt(9);
4739 brick_face_node_pt[10]=brick_el0_pt->node_pt(1);
4740 brick_face_node_pt[11]=brick_el1_pt->node_pt(3);
4742 brick_face_node_pt[12]=brick_el1_pt->node_pt(24);
4744 brick_face_node_pt[13]=brick_el0_pt->node_pt(10);
4745 brick_face_node_pt[14]=brick_el1_pt->node_pt(21);
4747 brick_face_node_pt[15]=brick_el1_pt->node_pt(12);
4748 brick_face_node_pt[16]=brick_el3_pt->node_pt(7);
4750 brick_face_node_pt[17]=brick_el3_pt->node_pt(4);
4751 brick_face_node_pt[18]=brick_el1_pt->node_pt(15);
4758 Vector<unsigned> translate(19);
4761 for (
unsigned i=0;i<19;i++)
4767 for (
unsigned j=0;j<19;j++)
4770 Node* brick_node_pt=brick_face_node_pt[translate[j]];
4773 Vector<double> s=s_face[j];
4774 Vector<double> zeta(2);
4775 Vector<double> x(3);
4776 face_el_pt->interpolated_zeta(s,zeta);
4777 face_el_pt->interpolated_x(s,x);
4781 double dist=sqrt(pow(brick_node_pt->x(0)-x[0],2)+
4782 pow(brick_node_pt->x(1)-x[1],2)+
4783 pow(brick_node_pt->x(2)-x[2],2));
4784 if (dist>BrickFromTetMeshHelper::Face_position_tolerance)
4786 std::ofstream brick0;
4787 std::ofstream brick1;
4788 std::ofstream brick2;
4789 std::ofstream brick3;
4790 brick0.open(
"full_brick0.dat");
4791 brick1.open(
"full_brick1.dat");
4792 brick2.open(
"full_brick2.dat");
4793 brick3.open(
"full_brick3.dat");
4794 for (
unsigned j=0;j<27;j++)
4796 brick0 << brick_el0_pt->node_pt(j)->x(0) <<
" " 4797 << brick_el0_pt->node_pt(j)->x(1) <<
" " 4798 << brick_el0_pt->node_pt(j)->x(2) <<
"\n";
4800 brick1 << brick_el1_pt->node_pt(j)->x(0) <<
" " 4801 << brick_el1_pt->node_pt(j)->x(1) <<
" " 4802 << brick_el1_pt->node_pt(j)->x(2) <<
"\n";
4804 brick2 << brick_el2_pt->node_pt(j)->x(0) <<
" " 4805 << brick_el2_pt->node_pt(j)->x(1) <<
" " 4806 << brick_el2_pt->node_pt(j)->x(2) <<
"\n";
4808 brick3 << brick_el3_pt->node_pt(j)->x(0) <<
" " 4809 << brick_el3_pt->node_pt(j)->x(1) <<
" " 4810 << brick_el3_pt->node_pt(j)->x(2) <<
"\n";
4817 std::ofstream full_face;
4818 full_face.open(
"full_face.dat");
4819 for (
unsigned j=0;j<6;j++)
4821 full_face << face_el_pt->node_pt(j)->x(0) <<
" " 4822 << face_el_pt->node_pt(j)->x(1) <<
" " 4823 << face_el_pt->node_pt(j)->x(2) <<
"\n";
4828 int normal_sign=face_el_pt->normal_sign();
4830 std::ostringstream error_stream;
4832 <<
"During assignment of boundary cordinates, the distance\n" 4833 <<
"between brick node and reference point in \n" 4834 <<
"triangular FaceElement is " << dist <<
" which \n" 4835 <<
"is bigger than the tolerance defined in \n" 4836 <<
"BrickFromTetMeshHelper::Face_position_tolerance=" 4837 << BrickFromTetMeshHelper::Face_position_tolerance <<
".\n" 4838 <<
"If this is tolerable, increase the tolerance \n" 4839 <<
"(it's defined in a namespace and therefore publically\n" 4840 <<
"accessible). If not, the Face may be inverted in which \n" 4841 <<
"case you should re-implement the translation scheme,\n" 4842 <<
"following the pattern used in the ThinLayerBrickOnTetMesh." 4843 <<
"\nThe required code fragements are already here but \n" 4844 <<
"the translation is the unit map.\n" 4845 <<
"To aid the diagnostics, the files full_brick[0-3].dat\n" 4846 <<
"contain the coordinates of the 27 nodes in the four\n" 4847 <<
"bricks associated with the current tet and full_face.dat\n" 4848 <<
"contains the coordinates of the 6 nodes in the FaceElement" 4849 <<
"\nfrom which the boundary coordinates are extracted.\n" 4850 <<
"FYI: The normal_sign of the face is: " << normal_sign
4852 throw OomphLibError(
4854 OOMPH_CURRENT_FUNCTION,
4855 OOMPH_EXCEPTION_LOCATION);
4860 add_boundary_node(b,brick_node_pt);
4861 brick_node_pt->set_coordinates_on_boundary(b,zeta);
4868 Boundary_element_pt[b].push_back(brick_el1_pt);
4869 Face_index_at_boundary[b].push_back(-2);
4870 Boundary_element_pt[b].push_back(brick_el2_pt);
4871 Face_index_at_boundary[b].push_back(-1);
4872 Boundary_element_pt[b].push_back(brick_el3_pt);
4873 Face_index_at_boundary[b].push_back(-2);
4877 Boundary_element_pt[b].push_back(brick_el0_pt);
4878 Face_index_at_boundary[b].push_back(-1);
4879 Boundary_element_pt[b].push_back(brick_el2_pt);
4880 Face_index_at_boundary[b].push_back(-2);
4881 Boundary_element_pt[b].push_back(brick_el3_pt);
4882 Face_index_at_boundary[b].push_back(-1);
4886 Boundary_element_pt[b].push_back(brick_el0_pt);
4887 Face_index_at_boundary[b].push_back(-3);
4888 Boundary_element_pt[b].push_back(brick_el1_pt);
4889 Face_index_at_boundary[b].push_back(-3);
4890 Boundary_element_pt[b].push_back(brick_el2_pt);
4891 Face_index_at_boundary[b].push_back(-3);
4895 Boundary_element_pt[b].push_back(brick_el0_pt);
4896 Face_index_at_boundary[b].push_back(-2);
4897 Boundary_element_pt[b].push_back(brick_el1_pt);
4898 Face_index_at_boundary[b].push_back(-1);
4899 Boundary_element_pt[b].push_back(brick_el3_pt);
4900 Face_index_at_boundary[b].push_back(-3);
4913 Lookup_for_elements_next_boundary_is_setup=
true;
4917 unsigned n_xda_boundaries=tet_mesh_pt->nxda_boundary();
4920 Boundary_id.resize(n_xda_boundaries);
4921 for (
unsigned xda_b=0;xda_b<n_xda_boundaries;xda_b++)
4923 Boundary_id[xda_b]=tet_mesh_pt->oomph_lib_boundary_ids(xda_b);
4928 for (
unsigned e=0;e<4;e++)
4930 for (
unsigned j=0;j<8;j++)
4932 delete dummy_q_el_pt[e]->node_pt(j);
4934 delete dummy_q_el_pt[e];
4943 template<
class ELEMENT>
4946 TimeStepper* time_stepper_pt)
4949 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3, 3);
4952 bool tet_mesh_is_solid_mesh=
false;
4953 if (dynamic_cast<SolidFiniteElement*>(tet_mesh_pt->element_pt(0))!=0)
4955 tet_mesh_is_solid_mesh=
true;
4962 Vector<Vector<double> > s_face(19);
4963 for (
unsigned i=0;i<19;i++)
4965 s_face[i].resize(2);
5040 s_face[i][0]=1.0/3.0;
5041 s_face[i][1]=1.0/3.0;
5048 s_face[i][0]=5.0/24.0;
5049 s_face[i][1]=5.0/24.0;
5053 s_face[i][0]=5.0/12.0;
5054 s_face[i][1]=5.0/12.0;
5060 s_face[i][1]=5.0/24.0;
5061 s_face[i][0]=7.0/12.0;
5065 s_face[i][1]=5.0/12.0;
5066 s_face[i][0]=1.0/6.0;
5073 s_face[i][0]=5.0/24.0;
5074 s_face[i][1]=7.0/12.0;
5078 s_face[i][0]=5.0/12.0;
5079 s_face[i][1]=1.0/6.0;
5086 unsigned nb=tet_mesh_pt->nboundary();
5090 Boundary_element_pt.resize(nb);
5091 Face_index_at_boundary.resize(nb);
5096 std::map<Node*,Node*> tet_node_node_pt;
5099 std::map<Edge,Node*> brick_edge_node_pt;
5102 std::map<TFace,Node*> tet_face_node_pt;
5106 Vector<DummyBrickElement*> dummy_q_el_pt(4);
5107 for (
unsigned e=0;e<4;e++)
5109 dummy_q_el_pt[e]=
new DummyBrickElement;
5110 for (
unsigned j=0;j<8;j++)
5112 dummy_q_el_pt[e]->construct_node(j);
5117 unsigned n_el_tet=tet_mesh_pt->nelement();
5118 for (
unsigned e_tet=0; e_tet<n_el_tet; e_tet++)
5121 TElement<3,3>* tet_el_pt=
dynamic_cast<TElement<3,3>*
>(
5122 tet_mesh_pt->element_pt(e_tet));
5127 std::ostringstream error_stream;
5129 <<
"BrickFromTetMesh can only built from tet mesh containing\n" 5130 <<
"ten-noded tets.\n";
5131 throw OomphLibError(
5133 OOMPH_CURRENT_FUNCTION,
5134 OOMPH_EXCEPTION_LOCATION);
5139 Node* centroid_node_pt=0;
5142 Node* top_mid_face_node0_pt=0;
5143 Node* right_mid_face_node0_pt=0;
5144 Node* back_mid_face_node0_pt=0;
5146 Node* top_mid_face_node1_pt=0;
5147 Node* right_mid_face_node1_pt=0;
5149 Node* top_mid_face_node2_pt=0;
5152 FiniteElement* brick_el0_pt=0;
5153 FiniteElement* brick_el1_pt=0;
5154 FiniteElement* brick_el2_pt=0;
5155 FiniteElement* brick_el3_pt=0;
5163 for (
unsigned j=0;j<8;j++)
5165 Node* nod_pt=dummy_q_el_pt[0]->node_pt(j);
5166 Vector<double> s_tet(3);
5167 Vector<double> x_tet(3);
5171 tet_el_pt->local_coordinate_of_node(0,s_tet);
5172 nod_pt->set_value(0,s_tet[0]);
5173 nod_pt->set_value(1,s_tet[1]);
5174 nod_pt->set_value(2,s_tet[2]);
5175 tet_el_pt->interpolated_x(s_tet,x_tet);
5176 nod_pt->x(0)=x_tet[0];
5177 nod_pt->x(1)=x_tet[1];
5178 nod_pt->x(2)=x_tet[2];
5181 tet_el_pt->local_coordinate_of_node(4,s_tet);
5182 nod_pt->set_value(0,s_tet[0]);
5183 nod_pt->set_value(1,s_tet[1]);
5184 nod_pt->set_value(2,s_tet[2]);
5185 tet_el_pt->interpolated_x(s_tet,x_tet);
5186 nod_pt->x(0)=x_tet[0];
5187 nod_pt->x(1)=x_tet[1];
5188 nod_pt->x(2)=x_tet[2];
5191 tet_el_pt->local_coordinate_of_node(6,s_tet);
5192 nod_pt->set_value(0,s_tet[0]);
5193 nod_pt->set_value(1,s_tet[1]);
5194 nod_pt->set_value(2,s_tet[2]);
5195 tet_el_pt->interpolated_x(s_tet,x_tet);
5196 nod_pt->x(0)=x_tet[0];
5197 nod_pt->x(1)=x_tet[1];
5198 nod_pt->x(2)=x_tet[2];
5206 nod_pt->set_value(0,s_tet[0]);
5207 nod_pt->set_value(1,s_tet[1]);
5208 nod_pt->set_value(2,s_tet[2]);
5209 tet_el_pt->interpolated_x(s_tet,x_tet);
5210 nod_pt->x(0)=x_tet[0];
5211 nod_pt->x(1)=x_tet[1];
5212 nod_pt->x(2)=x_tet[2];
5215 tet_el_pt->local_coordinate_of_node(5,s_tet);
5216 nod_pt->set_value(0,s_tet[0]);
5217 nod_pt->set_value(1,s_tet[1]);
5218 nod_pt->set_value(2,s_tet[2]);
5219 tet_el_pt->interpolated_x(s_tet,x_tet);
5220 nod_pt->x(0)=x_tet[0];
5221 nod_pt->x(1)=x_tet[1];
5222 nod_pt->x(2)=x_tet[2];
5230 nod_pt->set_value(0,s_tet[0]);
5231 nod_pt->set_value(1,s_tet[1]);
5232 nod_pt->set_value(2,s_tet[2]);
5233 tet_el_pt->interpolated_x(s_tet,x_tet);
5234 nod_pt->x(0)=x_tet[0];
5235 nod_pt->x(1)=x_tet[1];
5236 nod_pt->x(2)=x_tet[2];
5244 nod_pt->set_value(0,s_tet[0]);
5245 nod_pt->set_value(1,s_tet[1]);
5246 nod_pt->set_value(2,s_tet[2]);
5247 tet_el_pt->interpolated_x(s_tet,x_tet);
5248 nod_pt->x(0)=x_tet[0];
5249 nod_pt->x(1)=x_tet[1];
5250 nod_pt->x(2)=x_tet[2];
5257 nod_pt->set_value(0,s_tet[0]);
5258 nod_pt->set_value(1,s_tet[1]);
5259 nod_pt->set_value(2,s_tet[2]);
5260 tet_el_pt->interpolated_x(s_tet,x_tet);
5261 nod_pt->x(0)=x_tet[0];
5262 nod_pt->x(1)=x_tet[1];
5263 nod_pt->x(2)=x_tet[2];
5270 FiniteElement* el_pt=
new ELEMENT;
5272 Element_pt.push_back(el_pt);
5274 TFace face0(tet_el_pt->node_pt(0),
5275 tet_el_pt->node_pt(1),
5276 tet_el_pt->node_pt(2));
5278 TFace face1(tet_el_pt->node_pt(0),
5279 tet_el_pt->node_pt(2),
5280 tet_el_pt->node_pt(3));
5282 TFace face2(tet_el_pt->node_pt(0),
5283 tet_el_pt->node_pt(1),
5284 tet_el_pt->node_pt(3));
5288 Vector<Vector<unsigned> > tet_edge_node(3);
5289 tet_edge_node[0].resize(2);
5290 tet_edge_node[0][0]=4;
5291 tet_edge_node[0][1]=1;
5292 tet_edge_node[1].resize(2);
5293 tet_edge_node[1][0]=6;
5294 tet_edge_node[1][1]=3;
5295 tet_edge_node[2].resize(2);
5296 tet_edge_node[2][0]=5;
5297 tet_edge_node[2][1]=2;
5300 unsigned central_tet_vertex=0;
5302 Node* tet_node_pt=0;
5303 Node* old_node_pt=0;
5310 tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
5311 old_node_pt=tet_node_node_pt[tet_node_pt];
5314 Node* new_node_pt=0;
5315 if (tet_node_pt->is_on_boundary())
5317 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5321 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5323 tet_node_node_pt[tet_node_pt]=new_node_pt;
5324 Node_pt.push_back(new_node_pt);
5325 Vector<double> s(3);
5326 Vector<double> s_tet(3);
5327 Vector<double> x_tet(3);
5328 el_pt->local_coordinate_of_node(j,s);
5329 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5330 tet_el_pt->interpolated_x(s_tet,x_tet);
5331 new_node_pt->x(0)=x_tet[0];
5332 new_node_pt->x(1)=x_tet[1];
5333 new_node_pt->x(2)=x_tet[2];
5338 el_pt->node_pt(j)=old_node_pt;
5348 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
5349 old_node_pt=tet_node_node_pt[tet_node_pt];
5352 Node* new_node_pt=0;
5353 if (tet_node_pt->is_on_boundary())
5355 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5359 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5361 tet_node_node_pt[tet_node_pt]=new_node_pt;
5362 Node_pt.push_back(new_node_pt);
5363 Vector<double> s(3);
5364 Vector<double> s_tet(3);
5365 Vector<double> x_tet(3);
5366 el_pt->local_coordinate_of_node(j,s);
5367 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5368 tet_el_pt->interpolated_x(s_tet,x_tet);
5369 new_node_pt->x(0)=x_tet[0];
5370 new_node_pt->x(1)=x_tet[1];
5371 new_node_pt->x(2)=x_tet[2];
5376 el_pt->node_pt(j)=old_node_pt;
5386 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
5387 old_node_pt=tet_node_node_pt[tet_node_pt];
5390 Node* new_node_pt=0;
5391 if (tet_node_pt->is_on_boundary())
5393 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5397 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5399 tet_node_node_pt[tet_node_pt]=new_node_pt;
5400 Node_pt.push_back(new_node_pt);
5401 Vector<double> s(3);
5402 Vector<double> s_tet(3);
5403 Vector<double> x_tet(3);
5404 el_pt->local_coordinate_of_node(j,s);
5405 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5406 tet_el_pt->interpolated_x(s_tet,x_tet);
5407 new_node_pt->x(0)=x_tet[0];
5408 new_node_pt->x(1)=x_tet[1];
5409 new_node_pt->x(2)=x_tet[2];
5414 el_pt->node_pt(j)=old_node_pt;
5424 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
5425 old_node_pt=tet_node_node_pt[tet_node_pt];
5428 Node* new_node_pt=0;
5429 if (tet_node_pt->is_on_boundary())
5431 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5435 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5437 tet_node_node_pt[tet_node_pt]=new_node_pt;
5438 Node_pt.push_back(new_node_pt);
5439 Vector<double> s(3);
5440 Vector<double> s_tet(3);
5441 Vector<double> x_tet(3);
5442 el_pt->local_coordinate_of_node(j,s);
5443 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5444 tet_el_pt->interpolated_x(s_tet,x_tet);
5445 new_node_pt->x(0)=x_tet[0];
5446 new_node_pt->x(1)=x_tet[1];
5447 new_node_pt->x(2)=x_tet[2];
5452 el_pt->node_pt(j)=old_node_pt;
5464 old_node_pt=tet_face_node_pt[face0];
5467 Node* new_node_pt=0;
5468 if (face0.is_boundary_face())
5470 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5474 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5476 tet_face_node_pt[face0]=new_node_pt;
5477 Node_pt.push_back(new_node_pt);
5478 Vector<double> s(3);
5479 Vector<double> s_tet(3);
5480 Vector<double> x_tet(3);
5481 el_pt->local_coordinate_of_node(j,s);
5482 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5483 tet_el_pt->interpolated_x(s_tet,x_tet);
5484 new_node_pt->x(0)=x_tet[0];
5485 new_node_pt->x(1)=x_tet[1];
5486 new_node_pt->x(2)=x_tet[2];
5491 el_pt->node_pt(j)=old_node_pt;
5501 old_node_pt=tet_face_node_pt[face1];
5504 Node* new_node_pt=0;
5505 if (face1.is_boundary_face())
5507 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5511 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5513 tet_face_node_pt[face1]=new_node_pt;
5514 Node_pt.push_back(new_node_pt);
5515 Vector<double> s(3);
5516 Vector<double> s_tet(3);
5517 Vector<double> x_tet(3);
5518 el_pt->local_coordinate_of_node(j,s);
5519 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5520 tet_el_pt->interpolated_x(s_tet,x_tet);
5521 new_node_pt->x(0)=x_tet[0];
5522 new_node_pt->x(1)=x_tet[1];
5523 new_node_pt->x(2)=x_tet[2];
5528 el_pt->node_pt(j)=old_node_pt;
5538 old_node_pt=tet_face_node_pt[face2];
5541 Node* new_node_pt=0;
5542 if (face2.is_boundary_face())
5544 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5548 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5550 tet_face_node_pt[face2]=new_node_pt;
5551 Node_pt.push_back(new_node_pt);
5552 Vector<double> s(3);
5553 Vector<double> s_tet(3);
5554 Vector<double> x_tet(3);
5555 el_pt->local_coordinate_of_node(j,s);
5556 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5557 tet_el_pt->interpolated_x(s_tet,x_tet);
5558 new_node_pt->x(0)=x_tet[0];
5559 new_node_pt->x(1)=x_tet[1];
5560 new_node_pt->x(2)=x_tet[2];
5565 el_pt->node_pt(j)=old_node_pt;
5576 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5577 centroid_node_pt=new_node_pt;
5578 Node_pt.push_back(new_node_pt);
5579 Vector<double> s(3);
5580 Vector<double> s_tet(3);
5581 Vector<double> x_tet(3);
5582 el_pt->local_coordinate_of_node(j,s);
5583 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5584 tet_el_pt->interpolated_x(s_tet,x_tet);
5585 new_node_pt->x(0)=x_tet[0];
5586 new_node_pt->x(1)=x_tet[1];
5587 new_node_pt->x(2)=x_tet[2];
5598 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5599 Node_pt.push_back(new_node_pt);
5600 Vector<double> s(3);
5601 Vector<double> s_tet(3);
5602 Vector<double> x_tet(3);
5603 el_pt->local_coordinate_of_node(j,s);
5604 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5605 tet_el_pt->interpolated_x(s_tet,x_tet);
5606 new_node_pt->x(0)=x_tet[0];
5607 new_node_pt->x(1)=x_tet[1];
5608 new_node_pt->x(2)=x_tet[2];
5617 Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
5618 old_node_pt=brick_edge_node_pt[edge];
5621 Node* new_node_pt=0;
5622 if (edge.is_boundary_edge())
5624 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5628 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5630 brick_edge_node_pt[edge]=new_node_pt;
5631 Node_pt.push_back(new_node_pt);
5632 Vector<double> s(3);
5633 Vector<double> s_tet(3);
5634 Vector<double> x_tet(3);
5635 el_pt->local_coordinate_of_node(j,s);
5636 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5637 tet_el_pt->interpolated_x(s_tet,x_tet);
5638 new_node_pt->x(0)=x_tet[0];
5639 new_node_pt->x(1)=x_tet[1];
5640 new_node_pt->x(2)=x_tet[2];
5645 el_pt->node_pt(j)=old_node_pt;
5655 Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
5656 old_node_pt=brick_edge_node_pt[edge];
5659 Node* new_node_pt=0;
5660 if (edge.is_boundary_edge())
5662 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5666 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5668 brick_edge_node_pt[edge]=new_node_pt;
5669 Node_pt.push_back(new_node_pt);
5670 Vector<double> s(3);
5671 Vector<double> s_tet(3);
5672 Vector<double> x_tet(3);
5673 el_pt->local_coordinate_of_node(j,s);
5674 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5675 tet_el_pt->interpolated_x(s_tet,x_tet);
5676 new_node_pt->x(0)=x_tet[0];
5677 new_node_pt->x(1)=x_tet[1];
5678 new_node_pt->x(2)=x_tet[2];
5683 el_pt->node_pt(j)=old_node_pt;
5692 Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
5693 old_node_pt=brick_edge_node_pt[edge];
5696 Node* new_node_pt=0;
5697 if (edge.is_boundary_edge())
5699 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5703 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5705 brick_edge_node_pt[edge]=new_node_pt;
5706 Node_pt.push_back(new_node_pt);
5707 Vector<double> s(3);
5708 Vector<double> s_tet(3);
5709 Vector<double> x_tet(3);
5710 el_pt->local_coordinate_of_node(j,s);
5711 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5712 tet_el_pt->interpolated_x(s_tet,x_tet);
5713 new_node_pt->x(0)=x_tet[0];
5714 new_node_pt->x(1)=x_tet[1];
5715 new_node_pt->x(2)=x_tet[2];
5720 el_pt->node_pt(j)=old_node_pt;
5729 Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
5730 old_node_pt=brick_edge_node_pt[edge];
5733 Node* new_node_pt=0;
5734 if (edge.is_boundary_edge())
5736 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5740 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5742 brick_edge_node_pt[edge]=new_node_pt;
5743 Node_pt.push_back(new_node_pt);
5744 Vector<double> s(3);
5745 Vector<double> s_tet(3);
5746 Vector<double> x_tet(3);
5747 el_pt->local_coordinate_of_node(j,s);
5748 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5749 tet_el_pt->interpolated_x(s_tet,x_tet);
5750 new_node_pt->x(0)=x_tet[0];
5751 new_node_pt->x(1)=x_tet[1];
5752 new_node_pt->x(2)=x_tet[2];
5757 el_pt->node_pt(j)=old_node_pt;
5766 Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
5767 old_node_pt=brick_edge_node_pt[edge];
5770 Node* new_node_pt=0;
5771 if (edge.is_boundary_edge())
5773 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5777 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5779 brick_edge_node_pt[edge]=new_node_pt;
5780 Node_pt.push_back(new_node_pt);
5781 Vector<double> s(3);
5782 Vector<double> s_tet(3);
5783 Vector<double> x_tet(3);
5784 el_pt->local_coordinate_of_node(j,s);
5785 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5786 tet_el_pt->interpolated_x(s_tet,x_tet);
5787 new_node_pt->x(0)=x_tet[0];
5788 new_node_pt->x(1)=x_tet[1];
5789 new_node_pt->x(2)=x_tet[2];
5794 el_pt->node_pt(j)=old_node_pt;
5804 Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
5805 old_node_pt=brick_edge_node_pt[edge];
5808 Node* new_node_pt=0;
5809 if (edge.is_boundary_edge())
5811 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5815 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5817 brick_edge_node_pt[edge]=new_node_pt;
5818 Node_pt.push_back(new_node_pt);
5819 Vector<double> s(3);
5820 Vector<double> s_tet(3);
5821 Vector<double> x_tet(3);
5822 el_pt->local_coordinate_of_node(j,s);
5823 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5824 tet_el_pt->interpolated_x(s_tet,x_tet);
5825 new_node_pt->x(0)=x_tet[0];
5826 new_node_pt->x(1)=x_tet[1];
5827 new_node_pt->x(2)=x_tet[2];
5832 el_pt->node_pt(j)=old_node_pt;
5841 Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
5842 old_node_pt=brick_edge_node_pt[edge];
5845 Node* new_node_pt=0;
5846 if (edge.is_boundary_edge())
5848 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5852 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5854 brick_edge_node_pt[edge]=new_node_pt;
5855 Node_pt.push_back(new_node_pt);
5856 Vector<double> s(3);
5857 Vector<double> s_tet(3);
5858 Vector<double> x_tet(3);
5859 el_pt->local_coordinate_of_node(j,s);
5860 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5861 tet_el_pt->interpolated_x(s_tet,x_tet);
5862 new_node_pt->x(0)=x_tet[0];
5863 new_node_pt->x(1)=x_tet[1];
5864 new_node_pt->x(2)=x_tet[2];
5869 el_pt->node_pt(j)=old_node_pt;
5879 Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
5880 old_node_pt=brick_edge_node_pt[edge];
5883 Node* new_node_pt=0;
5884 if (edge.is_boundary_edge())
5886 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5890 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5892 brick_edge_node_pt[edge]=new_node_pt;
5893 Node_pt.push_back(new_node_pt);
5894 Vector<double> s(3);
5895 Vector<double> s_tet(3);
5896 Vector<double> x_tet(3);
5897 el_pt->local_coordinate_of_node(j,s);
5898 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5899 tet_el_pt->interpolated_x(s_tet,x_tet);
5900 new_node_pt->x(0)=x_tet[0];
5901 new_node_pt->x(1)=x_tet[1];
5902 new_node_pt->x(2)=x_tet[2];
5907 el_pt->node_pt(j)=old_node_pt;
5916 Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
5917 old_node_pt=brick_edge_node_pt[edge];
5920 Node* new_node_pt=0;
5921 if (edge.is_boundary_edge())
5923 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5927 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5929 brick_edge_node_pt[edge]=new_node_pt;
5930 Node_pt.push_back(new_node_pt);
5931 Vector<double> s(3);
5932 Vector<double> s_tet(3);
5933 Vector<double> x_tet(3);
5934 el_pt->local_coordinate_of_node(j,s);
5935 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5936 tet_el_pt->interpolated_x(s_tet,x_tet);
5937 new_node_pt->x(0)=x_tet[0];
5938 new_node_pt->x(1)=x_tet[1];
5939 new_node_pt->x(2)=x_tet[2];
5944 el_pt->node_pt(j)=old_node_pt;
5954 Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
5955 old_node_pt=brick_edge_node_pt[edge];
5958 Node* new_node_pt=0;
5959 if (edge.is_boundary_edge())
5961 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
5965 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
5967 brick_edge_node_pt[edge]=new_node_pt;
5968 Node_pt.push_back(new_node_pt);
5969 Vector<double> s(3);
5970 Vector<double> s_tet(3);
5971 Vector<double> x_tet(3);
5972 el_pt->local_coordinate_of_node(j,s);
5973 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
5974 tet_el_pt->interpolated_x(s_tet,x_tet);
5975 new_node_pt->x(0)=x_tet[0];
5976 new_node_pt->x(1)=x_tet[1];
5977 new_node_pt->x(2)=x_tet[2];
5982 el_pt->node_pt(j)=old_node_pt;
5992 Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
5993 old_node_pt=brick_edge_node_pt[edge];
5996 Node* new_node_pt=0;
5997 if (edge.is_boundary_edge())
5999 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6003 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6005 brick_edge_node_pt[edge]=new_node_pt;
6006 Node_pt.push_back(new_node_pt);
6007 Vector<double> s(3);
6008 Vector<double> s_tet(3);
6009 Vector<double> x_tet(3);
6010 el_pt->local_coordinate_of_node(j,s);
6011 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6012 tet_el_pt->interpolated_x(s_tet,x_tet);
6013 new_node_pt->x(0)=x_tet[0];
6014 new_node_pt->x(1)=x_tet[1];
6015 new_node_pt->x(2)=x_tet[2];
6020 el_pt->node_pt(j)=old_node_pt;
6030 Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
6031 old_node_pt=brick_edge_node_pt[edge];
6034 Node* new_node_pt=0;
6035 if (edge.is_boundary_edge())
6037 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6041 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6043 brick_edge_node_pt[edge]=new_node_pt;
6044 Node_pt.push_back(new_node_pt);
6045 Vector<double> s(3);
6046 Vector<double> s_tet(3);
6047 Vector<double> x_tet(3);
6048 el_pt->local_coordinate_of_node(j,s);
6049 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6050 tet_el_pt->interpolated_x(s_tet,x_tet);
6051 new_node_pt->x(0)=x_tet[0];
6052 new_node_pt->x(1)=x_tet[1];
6053 new_node_pt->x(2)=x_tet[2];
6058 el_pt->node_pt(j)=old_node_pt;
6069 TFace face(tet_el_pt->node_pt( central_tet_vertex),
6070 tet_el_pt->node_pt(tet_edge_node[0][0]),
6071 tet_el_pt->node_pt(tet_edge_node[2][0]));
6073 old_node_pt=tet_face_node_pt[face];
6076 Node* new_node_pt=0;
6077 if (face.is_boundary_face())
6079 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6083 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6085 tet_face_node_pt[face]=new_node_pt;
6086 Node_pt.push_back(new_node_pt);
6087 Vector<double> s(3);
6088 Vector<double> s_tet(3);
6089 Vector<double> x_tet(3);
6090 el_pt->local_coordinate_of_node(j,s);
6091 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6092 tet_el_pt->interpolated_x(s_tet,x_tet);
6093 new_node_pt->x(0)=x_tet[0];
6094 new_node_pt->x(1)=x_tet[1];
6095 new_node_pt->x(2)=x_tet[2];
6100 el_pt->node_pt(j)=old_node_pt;
6112 TFace face(tet_el_pt->node_pt(central_tet_vertex),
6113 tet_el_pt->node_pt(tet_edge_node[1][0]),
6114 tet_el_pt->node_pt(tet_edge_node[2][0]));
6116 old_node_pt=tet_face_node_pt[face];
6119 Node* new_node_pt=0;
6120 if (face.is_boundary_face())
6122 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6126 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6128 tet_face_node_pt[face]=new_node_pt;
6129 Node_pt.push_back(new_node_pt);
6130 Vector<double> s(3);
6131 Vector<double> s_tet(3);
6132 Vector<double> x_tet(3);
6133 el_pt->local_coordinate_of_node(j,s);
6134 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6135 tet_el_pt->interpolated_x(s_tet,x_tet);
6136 new_node_pt->x(0)=x_tet[0];
6137 new_node_pt->x(1)=x_tet[1];
6138 new_node_pt->x(2)=x_tet[2];
6143 el_pt->node_pt(j)=old_node_pt;
6155 TFace face(tet_el_pt->node_pt(central_tet_vertex),
6156 tet_el_pt->node_pt(tet_edge_node[0][0]),
6157 tet_el_pt->node_pt(tet_edge_node[1][0]));
6159 old_node_pt=tet_face_node_pt[face];
6162 Node* new_node_pt=0;
6163 if (face.is_boundary_face())
6165 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6169 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6171 tet_face_node_pt[face]=new_node_pt;
6172 Node_pt.push_back(new_node_pt);
6173 Vector<double> s(3);
6174 Vector<double> s_tet(3);
6175 Vector<double> x_tet(3);
6176 el_pt->local_coordinate_of_node(j,s);
6177 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6178 tet_el_pt->interpolated_x(s_tet,x_tet);
6179 new_node_pt->x(0)=x_tet[0];
6180 new_node_pt->x(1)=x_tet[1];
6181 new_node_pt->x(2)=x_tet[2];
6186 el_pt->node_pt(j)=old_node_pt;
6194 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6195 Node_pt.push_back(new_node_pt);
6196 Vector<double> s(3);
6197 Vector<double> s_tet(3);
6198 Vector<double> x_tet(3);
6199 el_pt->local_coordinate_of_node(j,s);
6200 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6201 top_mid_face_node0_pt=new_node_pt;
6202 tet_el_pt->interpolated_x(s_tet,x_tet);
6203 new_node_pt->x(0)=x_tet[0];
6204 new_node_pt->x(1)=x_tet[1];
6205 new_node_pt->x(2)=x_tet[2];
6213 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6214 Node_pt.push_back(new_node_pt);
6215 Vector<double> s(3);
6216 Vector<double> s_tet(3);
6217 Vector<double> x_tet(3);
6218 el_pt->local_coordinate_of_node(j,s);
6219 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6220 right_mid_face_node0_pt=new_node_pt;
6221 tet_el_pt->interpolated_x(s_tet,x_tet);
6222 new_node_pt->x(0)=x_tet[0];
6223 new_node_pt->x(1)=x_tet[1];
6224 new_node_pt->x(2)=x_tet[2];
6231 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6232 Node_pt.push_back(new_node_pt);
6233 Vector<double> s(3);
6234 Vector<double> s_tet(3);
6235 Vector<double> x_tet(3);
6236 el_pt->local_coordinate_of_node(j,s);
6237 dummy_q_el_pt[0]->interpolated_s_tet(s,s_tet);
6238 back_mid_face_node0_pt=new_node_pt;
6239 tet_el_pt->interpolated_x(s_tet,x_tet);
6240 new_node_pt->x(0)=x_tet[0];
6241 new_node_pt->x(1)=x_tet[1];
6242 new_node_pt->x(2)=x_tet[2];
6252 for (
unsigned j=0;j<8;j++)
6254 Node* nod_pt=dummy_q_el_pt[1]->node_pt(j);
6255 Vector<double> s_tet(3);
6256 Vector<double> x_tet(3);
6260 tet_el_pt->local_coordinate_of_node(1,s_tet);
6261 nod_pt->set_value(0,s_tet[0]);
6262 nod_pt->set_value(1,s_tet[1]);
6263 nod_pt->set_value(2,s_tet[2]);
6264 tet_el_pt->interpolated_x(s_tet,x_tet);
6265 nod_pt->x(0)=x_tet[0];
6266 nod_pt->x(1)=x_tet[1];
6267 nod_pt->x(2)=x_tet[2];
6270 tet_el_pt->local_coordinate_of_node(9,s_tet);
6271 nod_pt->set_value(0,s_tet[0]);
6272 nod_pt->set_value(1,s_tet[1]);
6273 nod_pt->set_value(2,s_tet[2]);
6274 tet_el_pt->interpolated_x(s_tet,x_tet);
6275 nod_pt->x(0)=x_tet[0];
6276 nod_pt->x(1)=x_tet[1];
6277 nod_pt->x(2)=x_tet[2];
6280 tet_el_pt->local_coordinate_of_node(4,s_tet);
6281 nod_pt->set_value(0,s_tet[0]);
6282 nod_pt->set_value(1,s_tet[1]);
6283 nod_pt->set_value(2,s_tet[2]);
6284 tet_el_pt->interpolated_x(s_tet,x_tet);
6285 nod_pt->x(0)=x_tet[0];
6286 nod_pt->x(1)=x_tet[1];
6287 nod_pt->x(2)=x_tet[2];
6295 nod_pt->set_value(0,s_tet[0]);
6296 nod_pt->set_value(1,s_tet[1]);
6297 nod_pt->set_value(2,s_tet[2]);
6298 tet_el_pt->interpolated_x(s_tet,x_tet);
6299 nod_pt->x(0)=x_tet[0];
6300 nod_pt->x(1)=x_tet[1];
6301 nod_pt->x(2)=x_tet[2];
6304 tet_el_pt->local_coordinate_of_node(7,s_tet);
6305 nod_pt->set_value(0,s_tet[0]);
6306 nod_pt->set_value(1,s_tet[1]);
6307 nod_pt->set_value(2,s_tet[2]);
6308 tet_el_pt->interpolated_x(s_tet,x_tet);
6309 nod_pt->x(0)=x_tet[0];
6310 nod_pt->x(1)=x_tet[1];
6311 nod_pt->x(2)=x_tet[2];
6319 nod_pt->set_value(0,s_tet[0]);
6320 nod_pt->set_value(1,s_tet[1]);
6321 nod_pt->set_value(2,s_tet[2]);
6322 tet_el_pt->interpolated_x(s_tet,x_tet);
6323 nod_pt->x(0)=x_tet[0];
6324 nod_pt->x(1)=x_tet[1];
6325 nod_pt->x(2)=x_tet[2];
6333 nod_pt->set_value(0,s_tet[0]);
6334 nod_pt->set_value(1,s_tet[1]);
6335 nod_pt->set_value(2,s_tet[2]);
6336 tet_el_pt->interpolated_x(s_tet,x_tet);
6337 nod_pt->x(0)=x_tet[0];
6338 nod_pt->x(1)=x_tet[1];
6339 nod_pt->x(2)=x_tet[2];
6346 nod_pt->set_value(0,s_tet[0]);
6347 nod_pt->set_value(1,s_tet[1]);
6348 nod_pt->set_value(2,s_tet[2]);
6349 tet_el_pt->interpolated_x(s_tet,x_tet);
6350 nod_pt->x(0)=x_tet[0];
6351 nod_pt->x(1)=x_tet[1];
6352 nod_pt->x(2)=x_tet[2];
6359 FiniteElement* el_pt=
new ELEMENT;
6361 Element_pt.push_back(el_pt);
6363 TFace face0(tet_el_pt->node_pt(1),
6364 tet_el_pt->node_pt(3),
6365 tet_el_pt->node_pt(2));
6367 TFace face1(tet_el_pt->node_pt(1),
6368 tet_el_pt->node_pt(0),
6369 tet_el_pt->node_pt(2));
6371 TFace face2(tet_el_pt->node_pt(1),
6372 tet_el_pt->node_pt(0),
6373 tet_el_pt->node_pt(3));
6376 Vector<Vector<unsigned> > tet_edge_node(3);
6377 tet_edge_node[0].resize(2);
6378 tet_edge_node[0][0]=9;
6379 tet_edge_node[0][1]=3;
6380 tet_edge_node[1].resize(2);
6381 tet_edge_node[1][0]=4;
6382 tet_edge_node[1][1]=0;
6383 tet_edge_node[2].resize(2);
6384 tet_edge_node[2][0]=7;
6385 tet_edge_node[2][1]=2;
6388 unsigned central_tet_vertex=1;
6390 Node* tet_node_pt=0;
6391 Node* old_node_pt=0;
6398 tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
6399 old_node_pt=tet_node_node_pt[tet_node_pt];
6402 Node* new_node_pt=0;
6403 if (tet_node_pt->is_on_boundary())
6405 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6409 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6411 tet_node_node_pt[tet_node_pt]=new_node_pt;
6412 Node_pt.push_back(new_node_pt);
6413 Vector<double> s(3);
6414 Vector<double> s_tet(3);
6415 Vector<double> x_tet(3);
6416 el_pt->local_coordinate_of_node(j,s);
6417 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6418 tet_el_pt->interpolated_x(s_tet,x_tet);
6419 new_node_pt->x(0)=x_tet[0];
6420 new_node_pt->x(1)=x_tet[1];
6421 new_node_pt->x(2)=x_tet[2];
6426 el_pt->node_pt(j)=old_node_pt;
6436 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
6437 old_node_pt=tet_node_node_pt[tet_node_pt];
6440 Node* new_node_pt=0;
6441 if (tet_node_pt->is_on_boundary())
6443 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6447 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6449 tet_node_node_pt[tet_node_pt]=new_node_pt;
6450 Node_pt.push_back(new_node_pt);
6451 Vector<double> s(3);
6452 Vector<double> s_tet(3);
6453 Vector<double> x_tet(3);
6454 el_pt->local_coordinate_of_node(j,s);
6455 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6456 tet_el_pt->interpolated_x(s_tet,x_tet);
6457 new_node_pt->x(0)=x_tet[0];
6458 new_node_pt->x(1)=x_tet[1];
6459 new_node_pt->x(2)=x_tet[2];
6464 el_pt->node_pt(j)=old_node_pt;
6474 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
6475 old_node_pt=tet_node_node_pt[tet_node_pt];
6478 Node* new_node_pt=0;
6479 if (tet_node_pt->is_on_boundary())
6481 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6485 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6487 tet_node_node_pt[tet_node_pt]=new_node_pt;
6488 Node_pt.push_back(new_node_pt);
6489 Vector<double> s(3);
6490 Vector<double> s_tet(3);
6491 Vector<double> x_tet(3);
6492 el_pt->local_coordinate_of_node(j,s);
6493 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6494 tet_el_pt->interpolated_x(s_tet,x_tet);
6495 new_node_pt->x(0)=x_tet[0];
6496 new_node_pt->x(1)=x_tet[1];
6497 new_node_pt->x(2)=x_tet[2];
6502 el_pt->node_pt(j)=old_node_pt;
6512 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
6513 old_node_pt=tet_node_node_pt[tet_node_pt];
6516 Node* new_node_pt=0;
6517 if (tet_node_pt->is_on_boundary())
6519 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6523 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6525 tet_node_node_pt[tet_node_pt]=new_node_pt;
6526 Node_pt.push_back(new_node_pt);
6527 Vector<double> s(3);
6528 Vector<double> s_tet(3);
6529 Vector<double> x_tet(3);
6530 el_pt->local_coordinate_of_node(j,s);
6531 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6532 tet_el_pt->interpolated_x(s_tet,x_tet);
6533 new_node_pt->x(0)=x_tet[0];
6534 new_node_pt->x(1)=x_tet[1];
6535 new_node_pt->x(2)=x_tet[2];
6540 el_pt->node_pt(j)=old_node_pt;
6551 old_node_pt=tet_face_node_pt[face0];
6554 Node* new_node_pt=0;
6555 if (face0.is_boundary_face())
6557 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6561 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6563 tet_face_node_pt[face0]=new_node_pt;
6564 Node_pt.push_back(new_node_pt);
6565 Vector<double> s(3);
6566 Vector<double> s_tet(3);
6567 Vector<double> x_tet(3);
6568 el_pt->local_coordinate_of_node(j,s);
6569 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6570 tet_el_pt->interpolated_x(s_tet,x_tet);
6571 new_node_pt->x(0)=x_tet[0];
6572 new_node_pt->x(1)=x_tet[1];
6573 new_node_pt->x(2)=x_tet[2];
6578 el_pt->node_pt(j)=old_node_pt;
6587 old_node_pt=tet_face_node_pt[face1];
6590 Node* new_node_pt=0;
6591 if (face1.is_boundary_face())
6593 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6597 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6599 tet_face_node_pt[face1]=new_node_pt;
6600 Node_pt.push_back(new_node_pt);
6601 Vector<double> s(3);
6602 Vector<double> s_tet(3);
6603 Vector<double> x_tet(3);
6604 el_pt->local_coordinate_of_node(j,s);
6605 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6606 tet_el_pt->interpolated_x(s_tet,x_tet);
6607 new_node_pt->x(0)=x_tet[0];
6608 new_node_pt->x(1)=x_tet[1];
6609 new_node_pt->x(2)=x_tet[2];
6614 el_pt->node_pt(j)=old_node_pt;
6623 old_node_pt=tet_face_node_pt[face2];
6626 Node* new_node_pt=0;
6627 if (face2.is_boundary_face())
6629 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6633 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6635 tet_face_node_pt[face2]=new_node_pt;
6636 Node_pt.push_back(new_node_pt);
6637 Vector<double> s(3);
6638 Vector<double> s_tet(3);
6639 Vector<double> x_tet(3);
6640 el_pt->local_coordinate_of_node(j,s);
6641 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6642 tet_el_pt->interpolated_x(s_tet,x_tet);
6643 new_node_pt->x(0)=x_tet[0];
6644 new_node_pt->x(1)=x_tet[1];
6645 new_node_pt->x(2)=x_tet[2];
6650 el_pt->node_pt(j)=old_node_pt;
6661 el_pt->node_pt(j)=centroid_node_pt;
6671 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6672 Node_pt.push_back(new_node_pt);
6673 Vector<double> s(3);
6674 Vector<double> s_tet(3);
6675 Vector<double> x_tet(3);
6676 el_pt->local_coordinate_of_node(j,s);
6677 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6678 tet_el_pt->interpolated_x(s_tet,x_tet);
6679 new_node_pt->x(0)=x_tet[0];
6680 new_node_pt->x(1)=x_tet[1];
6681 new_node_pt->x(2)=x_tet[2];
6691 Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
6692 old_node_pt=brick_edge_node_pt[edge];
6695 Node* new_node_pt=0;
6696 if (edge.is_boundary_edge())
6698 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6702 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6704 brick_edge_node_pt[edge]=new_node_pt;
6705 Node_pt.push_back(new_node_pt);
6706 Vector<double> s(3);
6707 Vector<double> s_tet(3);
6708 Vector<double> x_tet(3);
6709 el_pt->local_coordinate_of_node(j,s);
6710 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6711 tet_el_pt->interpolated_x(s_tet,x_tet);
6712 new_node_pt->x(0)=x_tet[0];
6713 new_node_pt->x(1)=x_tet[1];
6714 new_node_pt->x(2)=x_tet[2];
6719 el_pt->node_pt(j)=old_node_pt;
6729 Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
6730 old_node_pt=brick_edge_node_pt[edge];
6733 Node* new_node_pt=0;
6734 if (edge.is_boundary_edge())
6736 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6740 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6742 brick_edge_node_pt[edge]=new_node_pt;
6743 Node_pt.push_back(new_node_pt);
6744 Vector<double> s(3);
6745 Vector<double> s_tet(3);
6746 Vector<double> x_tet(3);
6747 el_pt->local_coordinate_of_node(j,s);
6748 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6749 tet_el_pt->interpolated_x(s_tet,x_tet);
6750 new_node_pt->x(0)=x_tet[0];
6751 new_node_pt->x(1)=x_tet[1];
6752 new_node_pt->x(2)=x_tet[2];
6757 el_pt->node_pt(j)=old_node_pt;
6767 Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
6768 old_node_pt=brick_edge_node_pt[edge];
6771 Node* new_node_pt=0;
6772 if (edge.is_boundary_edge())
6774 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6778 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6780 brick_edge_node_pt[edge]=new_node_pt;
6781 Node_pt.push_back(new_node_pt);
6782 Vector<double> s(3);
6783 Vector<double> s_tet(3);
6784 Vector<double> x_tet(3);
6785 el_pt->local_coordinate_of_node(j,s);
6786 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6787 tet_el_pt->interpolated_x(s_tet,x_tet);
6788 new_node_pt->x(0)=x_tet[0];
6789 new_node_pt->x(1)=x_tet[1];
6790 new_node_pt->x(2)=x_tet[2];
6795 el_pt->node_pt(j)=old_node_pt;
6804 Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
6805 old_node_pt=brick_edge_node_pt[edge];
6808 Node* new_node_pt=0;
6809 if (edge.is_boundary_edge())
6811 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6815 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6817 brick_edge_node_pt[edge]=new_node_pt;
6818 Node_pt.push_back(new_node_pt);
6819 Vector<double> s(3);
6820 Vector<double> s_tet(3);
6821 Vector<double> x_tet(3);
6822 el_pt->local_coordinate_of_node(j,s);
6823 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6824 tet_el_pt->interpolated_x(s_tet,x_tet);
6825 new_node_pt->x(0)=x_tet[0];
6826 new_node_pt->x(1)=x_tet[1];
6827 new_node_pt->x(2)=x_tet[2];
6832 el_pt->node_pt(j)=old_node_pt;
6841 Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
6842 old_node_pt=brick_edge_node_pt[edge];
6845 Node* new_node_pt=0;
6846 if (edge.is_boundary_edge())
6848 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6852 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6854 brick_edge_node_pt[edge]=new_node_pt;
6855 Node_pt.push_back(new_node_pt);
6856 Vector<double> s(3);
6857 Vector<double> s_tet(3);
6858 Vector<double> x_tet(3);
6859 el_pt->local_coordinate_of_node(j,s);
6860 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6861 tet_el_pt->interpolated_x(s_tet,x_tet);
6862 new_node_pt->x(0)=x_tet[0];
6863 new_node_pt->x(1)=x_tet[1];
6864 new_node_pt->x(2)=x_tet[2];
6869 el_pt->node_pt(j)=old_node_pt;
6879 Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
6880 old_node_pt=brick_edge_node_pt[edge];
6883 Node* new_node_pt=0;
6884 if (edge.is_boundary_edge())
6886 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6890 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6892 brick_edge_node_pt[edge]=new_node_pt;
6893 Node_pt.push_back(new_node_pt);
6894 Vector<double> s(3);
6895 Vector<double> s_tet(3);
6896 Vector<double> x_tet(3);
6897 el_pt->local_coordinate_of_node(j,s);
6898 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6899 tet_el_pt->interpolated_x(s_tet,x_tet);
6900 new_node_pt->x(0)=x_tet[0];
6901 new_node_pt->x(1)=x_tet[1];
6902 new_node_pt->x(2)=x_tet[2];
6907 el_pt->node_pt(j)=old_node_pt;
6916 Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
6917 old_node_pt=brick_edge_node_pt[edge];
6920 Node* new_node_pt=0;
6921 if (edge.is_boundary_edge())
6923 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6927 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6929 brick_edge_node_pt[edge]=new_node_pt;
6930 Node_pt.push_back(new_node_pt);
6931 Vector<double> s(3);
6932 Vector<double> s_tet(3);
6933 Vector<double> x_tet(3);
6934 el_pt->local_coordinate_of_node(j,s);
6935 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6936 tet_el_pt->interpolated_x(s_tet,x_tet);
6937 new_node_pt->x(0)=x_tet[0];
6938 new_node_pt->x(1)=x_tet[1];
6939 new_node_pt->x(2)=x_tet[2];
6944 el_pt->node_pt(j)=old_node_pt;
6954 Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
6955 old_node_pt=brick_edge_node_pt[edge];
6958 Node* new_node_pt=0;
6959 if (edge.is_boundary_edge())
6961 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
6965 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
6967 brick_edge_node_pt[edge]=new_node_pt;
6968 Node_pt.push_back(new_node_pt);
6969 Vector<double> s(3);
6970 Vector<double> s_tet(3);
6971 Vector<double> x_tet(3);
6972 el_pt->local_coordinate_of_node(j,s);
6973 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
6974 tet_el_pt->interpolated_x(s_tet,x_tet);
6975 new_node_pt->x(0)=x_tet[0];
6976 new_node_pt->x(1)=x_tet[1];
6977 new_node_pt->x(2)=x_tet[2];
6982 el_pt->node_pt(j)=old_node_pt;
6991 Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
6992 old_node_pt=brick_edge_node_pt[edge];
6995 Node* new_node_pt=0;
6996 if (edge.is_boundary_edge())
6998 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7002 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7004 brick_edge_node_pt[edge]=new_node_pt;
7005 Node_pt.push_back(new_node_pt);
7006 Vector<double> s(3);
7007 Vector<double> s_tet(3);
7008 Vector<double> x_tet(3);
7009 el_pt->local_coordinate_of_node(j,s);
7010 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7011 tet_el_pt->interpolated_x(s_tet,x_tet);
7012 new_node_pt->x(0)=x_tet[0];
7013 new_node_pt->x(1)=x_tet[1];
7014 new_node_pt->x(2)=x_tet[2];
7019 el_pt->node_pt(j)=old_node_pt;
7029 Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
7030 old_node_pt=brick_edge_node_pt[edge];
7033 Node* new_node_pt=0;
7034 if (edge.is_boundary_edge())
7036 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7040 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7042 brick_edge_node_pt[edge]=new_node_pt;
7043 Node_pt.push_back(new_node_pt);
7044 Vector<double> s(3);
7045 Vector<double> s_tet(3);
7046 Vector<double> x_tet(3);
7047 el_pt->local_coordinate_of_node(j,s);
7048 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7049 tet_el_pt->interpolated_x(s_tet,x_tet);
7050 new_node_pt->x(0)=x_tet[0];
7051 new_node_pt->x(1)=x_tet[1];
7052 new_node_pt->x(2)=x_tet[2];
7057 el_pt->node_pt(j)=old_node_pt;
7067 Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
7068 old_node_pt=brick_edge_node_pt[edge];
7071 Node* new_node_pt=0;
7072 if (edge.is_boundary_edge())
7074 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7078 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7080 brick_edge_node_pt[edge]=new_node_pt;
7081 Node_pt.push_back(new_node_pt);
7082 Vector<double> s(3);
7083 Vector<double> s_tet(3);
7084 Vector<double> x_tet(3);
7085 el_pt->local_coordinate_of_node(j,s);
7086 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7087 tet_el_pt->interpolated_x(s_tet,x_tet);
7088 new_node_pt->x(0)=x_tet[0];
7089 new_node_pt->x(1)=x_tet[1];
7090 new_node_pt->x(2)=x_tet[2];
7095 el_pt->node_pt(j)=old_node_pt;
7105 Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
7106 old_node_pt=brick_edge_node_pt[edge];
7109 Node* new_node_pt=0;
7110 if (edge.is_boundary_edge())
7112 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7116 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7118 brick_edge_node_pt[edge]=new_node_pt;
7119 Node_pt.push_back(new_node_pt);
7120 Vector<double> s(3);
7121 Vector<double> s_tet(3);
7122 Vector<double> x_tet(3);
7123 el_pt->local_coordinate_of_node(j,s);
7124 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7125 tet_el_pt->interpolated_x(s_tet,x_tet);
7126 new_node_pt->x(0)=x_tet[0];
7127 new_node_pt->x(1)=x_tet[1];
7128 new_node_pt->x(2)=x_tet[2];
7133 el_pt->node_pt(j)=old_node_pt;
7144 TFace face(tet_el_pt->node_pt( central_tet_vertex),
7145 tet_el_pt->node_pt(tet_edge_node[0][0]),
7146 tet_el_pt->node_pt(tet_edge_node[2][0]));
7148 old_node_pt=tet_face_node_pt[face];
7151 Node* new_node_pt=0;
7152 if (face.is_boundary_face())
7154 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7158 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7160 tet_face_node_pt[face]=new_node_pt;
7161 Node_pt.push_back(new_node_pt);
7162 Vector<double> s(3);
7163 Vector<double> s_tet(3);
7164 Vector<double> x_tet(3);
7165 el_pt->local_coordinate_of_node(j,s);
7166 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7167 tet_el_pt->interpolated_x(s_tet,x_tet);
7168 new_node_pt->x(0)=x_tet[0];
7169 new_node_pt->x(1)=x_tet[1];
7170 new_node_pt->x(2)=x_tet[2];
7175 el_pt->node_pt(j)=old_node_pt;
7187 TFace face(tet_el_pt->node_pt(central_tet_vertex),
7188 tet_el_pt->node_pt(tet_edge_node[1][0]),
7189 tet_el_pt->node_pt(tet_edge_node[2][0]));
7191 old_node_pt=tet_face_node_pt[face];
7194 Node* new_node_pt=0;
7195 if (face.is_boundary_face())
7197 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7201 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7203 tet_face_node_pt[face]=new_node_pt;
7204 Node_pt.push_back(new_node_pt);
7205 Vector<double> s(3);
7206 Vector<double> s_tet(3);
7207 Vector<double> x_tet(3);
7208 el_pt->local_coordinate_of_node(j,s);
7209 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7210 tet_el_pt->interpolated_x(s_tet,x_tet);
7211 new_node_pt->x(0)=x_tet[0];
7212 new_node_pt->x(1)=x_tet[1];
7213 new_node_pt->x(2)=x_tet[2];
7218 el_pt->node_pt(j)=old_node_pt;
7230 TFace face(tet_el_pt->node_pt(central_tet_vertex),
7231 tet_el_pt->node_pt(tet_edge_node[0][0]),
7232 tet_el_pt->node_pt(tet_edge_node[1][0]));
7234 old_node_pt=tet_face_node_pt[face];
7237 Node* new_node_pt=0;
7238 if (face.is_boundary_face())
7240 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7244 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7246 tet_face_node_pt[face]=new_node_pt;
7247 Node_pt.push_back(new_node_pt);
7248 Vector<double> s(3);
7249 Vector<double> s_tet(3);
7250 Vector<double> x_tet(3);
7251 el_pt->local_coordinate_of_node(j,s);
7252 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7253 tet_el_pt->interpolated_x(s_tet,x_tet);
7254 new_node_pt->x(0)=x_tet[0];
7255 new_node_pt->x(1)=x_tet[1];
7256 new_node_pt->x(2)=x_tet[2];
7261 el_pt->node_pt(j)=old_node_pt;
7269 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7270 Node_pt.push_back(new_node_pt);
7271 Vector<double> s(3);
7272 Vector<double> s_tet(3);
7273 Vector<double> x_tet(3);
7274 el_pt->local_coordinate_of_node(j,s);
7275 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7276 top_mid_face_node1_pt=new_node_pt;
7277 tet_el_pt->interpolated_x(s_tet,x_tet);
7278 new_node_pt->x(0)=x_tet[0];
7279 new_node_pt->x(1)=x_tet[1];
7280 new_node_pt->x(2)=x_tet[2];
7288 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7289 Node_pt.push_back(new_node_pt);
7290 Vector<double> s(3);
7291 Vector<double> s_tet(3);
7292 Vector<double> x_tet(3);
7293 el_pt->local_coordinate_of_node(j,s);
7294 dummy_q_el_pt[1]->interpolated_s_tet(s,s_tet);
7295 right_mid_face_node1_pt=new_node_pt;
7296 tet_el_pt->interpolated_x(s_tet,x_tet);
7297 new_node_pt->x(0)=x_tet[0];
7298 new_node_pt->x(1)=x_tet[1];
7299 new_node_pt->x(2)=x_tet[2];
7308 el_pt->node_pt(j)=right_mid_face_node0_pt;
7317 for (
unsigned j=0;j<8;j++)
7319 Node* nod_pt=dummy_q_el_pt[2]->node_pt(j);
7320 Vector<double> s_tet(3);
7321 Vector<double> x_tet(3);
7325 tet_el_pt->local_coordinate_of_node(3,s_tet);
7326 nod_pt->set_value(0,s_tet[0]);
7327 nod_pt->set_value(1,s_tet[1]);
7328 nod_pt->set_value(2,s_tet[2]);
7329 tet_el_pt->interpolated_x(s_tet,x_tet);
7330 nod_pt->x(0)=x_tet[0];
7331 nod_pt->x(1)=x_tet[1];
7332 nod_pt->x(2)=x_tet[2];
7335 tet_el_pt->local_coordinate_of_node(6,s_tet);
7336 nod_pt->set_value(0,s_tet[0]);
7337 nod_pt->set_value(1,s_tet[1]);
7338 nod_pt->set_value(2,s_tet[2]);
7339 tet_el_pt->interpolated_x(s_tet,x_tet);
7340 nod_pt->x(0)=x_tet[0];
7341 nod_pt->x(1)=x_tet[1];
7342 nod_pt->x(2)=x_tet[2];
7345 tet_el_pt->local_coordinate_of_node(9,s_tet);
7346 nod_pt->set_value(0,s_tet[0]);
7347 nod_pt->set_value(1,s_tet[1]);
7348 nod_pt->set_value(2,s_tet[2]);
7349 tet_el_pt->interpolated_x(s_tet,x_tet);
7350 nod_pt->x(0)=x_tet[0];
7351 nod_pt->x(1)=x_tet[1];
7352 nod_pt->x(2)=x_tet[2];
7360 nod_pt->set_value(0,s_tet[0]);
7361 nod_pt->set_value(1,s_tet[1]);
7362 nod_pt->set_value(2,s_tet[2]);
7363 tet_el_pt->interpolated_x(s_tet,x_tet);
7364 nod_pt->x(0)=x_tet[0];
7365 nod_pt->x(1)=x_tet[1];
7366 nod_pt->x(2)=x_tet[2];
7369 tet_el_pt->local_coordinate_of_node(8,s_tet);
7370 nod_pt->set_value(0,s_tet[0]);
7371 nod_pt->set_value(1,s_tet[1]);
7372 nod_pt->set_value(2,s_tet[2]);
7373 tet_el_pt->interpolated_x(s_tet,x_tet);
7374 nod_pt->x(0)=x_tet[0];
7375 nod_pt->x(1)=x_tet[1];
7376 nod_pt->x(2)=x_tet[2];
7384 nod_pt->set_value(0,s_tet[0]);
7385 nod_pt->set_value(1,s_tet[1]);
7386 nod_pt->set_value(2,s_tet[2]);
7387 tet_el_pt->interpolated_x(s_tet,x_tet);
7388 nod_pt->x(0)=x_tet[0];
7389 nod_pt->x(1)=x_tet[1];
7390 nod_pt->x(2)=x_tet[2];
7398 nod_pt->set_value(0,s_tet[0]);
7399 nod_pt->set_value(1,s_tet[1]);
7400 nod_pt->set_value(2,s_tet[2]);
7401 tet_el_pt->interpolated_x(s_tet,x_tet);
7402 nod_pt->x(0)=x_tet[0];
7403 nod_pt->x(1)=x_tet[1];
7404 nod_pt->x(2)=x_tet[2];
7411 nod_pt->set_value(0,s_tet[0]);
7412 nod_pt->set_value(1,s_tet[1]);
7413 nod_pt->set_value(2,s_tet[2]);
7414 tet_el_pt->interpolated_x(s_tet,x_tet);
7415 nod_pt->x(0)=x_tet[0];
7416 nod_pt->x(1)=x_tet[1];
7417 nod_pt->x(2)=x_tet[2];
7424 FiniteElement* el_pt=
new ELEMENT;
7426 Element_pt.push_back(el_pt);
7428 TFace face0(tet_el_pt->node_pt(3),
7429 tet_el_pt->node_pt(0),
7430 tet_el_pt->node_pt(2));
7432 TFace face1(tet_el_pt->node_pt(3),
7433 tet_el_pt->node_pt(1),
7434 tet_el_pt->node_pt(2));
7436 TFace face2(tet_el_pt->node_pt(3),
7437 tet_el_pt->node_pt(1),
7438 tet_el_pt->node_pt(0));
7441 Vector<Vector<unsigned> > tet_edge_node(3);
7442 tet_edge_node[0].resize(2);
7443 tet_edge_node[0][0]=6;
7444 tet_edge_node[0][1]=0;
7445 tet_edge_node[1].resize(2);
7446 tet_edge_node[1][0]=9;
7447 tet_edge_node[1][1]=1;
7448 tet_edge_node[2].resize(2);
7449 tet_edge_node[2][0]=8;
7450 tet_edge_node[2][1]=2;
7453 unsigned central_tet_vertex=3;
7455 Node* tet_node_pt=0;
7456 Node* old_node_pt=0;
7463 tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
7464 old_node_pt=tet_node_node_pt[tet_node_pt];
7467 Node* new_node_pt=0;
7468 if (tet_node_pt->is_on_boundary())
7470 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7474 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7476 tet_node_node_pt[tet_node_pt]=new_node_pt;
7477 Node_pt.push_back(new_node_pt);
7478 Vector<double> s(3);
7479 Vector<double> s_tet(3);
7480 Vector<double> x_tet(3);
7481 el_pt->local_coordinate_of_node(j,s);
7482 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7483 tet_el_pt->interpolated_x(s_tet,x_tet);
7484 new_node_pt->x(0)=x_tet[0];
7485 new_node_pt->x(1)=x_tet[1];
7486 new_node_pt->x(2)=x_tet[2];
7491 el_pt->node_pt(j)=old_node_pt;
7501 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
7502 old_node_pt=tet_node_node_pt[tet_node_pt];
7505 Node* new_node_pt=0;
7506 if (tet_node_pt->is_on_boundary())
7508 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7512 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7514 tet_node_node_pt[tet_node_pt]=new_node_pt;
7515 Node_pt.push_back(new_node_pt);
7516 Vector<double> s(3);
7517 Vector<double> s_tet(3);
7518 Vector<double> x_tet(3);
7519 el_pt->local_coordinate_of_node(j,s);
7520 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7521 tet_el_pt->interpolated_x(s_tet,x_tet);
7522 new_node_pt->x(0)=x_tet[0];
7523 new_node_pt->x(1)=x_tet[1];
7524 new_node_pt->x(2)=x_tet[2];
7529 el_pt->node_pt(j)=old_node_pt;
7539 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
7540 old_node_pt=tet_node_node_pt[tet_node_pt];
7543 Node* new_node_pt=0;
7544 if (tet_node_pt->is_on_boundary())
7546 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7550 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7552 tet_node_node_pt[tet_node_pt]=new_node_pt;
7553 Node_pt.push_back(new_node_pt);
7554 Vector<double> s(3);
7555 Vector<double> s_tet(3);
7556 Vector<double> x_tet(3);
7557 el_pt->local_coordinate_of_node(j,s);
7558 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7559 tet_el_pt->interpolated_x(s_tet,x_tet);
7560 new_node_pt->x(0)=x_tet[0];
7561 new_node_pt->x(1)=x_tet[1];
7562 new_node_pt->x(2)=x_tet[2];
7567 el_pt->node_pt(j)=old_node_pt;
7577 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
7578 old_node_pt=tet_node_node_pt[tet_node_pt];
7581 Node* new_node_pt=0;
7582 if (tet_node_pt->is_on_boundary())
7584 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7588 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7590 tet_node_node_pt[tet_node_pt]=new_node_pt;
7591 Node_pt.push_back(new_node_pt);
7592 Vector<double> s(3);
7593 Vector<double> s_tet(3);
7594 Vector<double> x_tet(3);
7595 el_pt->local_coordinate_of_node(j,s);
7596 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7597 tet_el_pt->interpolated_x(s_tet,x_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];
7605 el_pt->node_pt(j)=old_node_pt;
7616 old_node_pt=tet_face_node_pt[face0];
7619 Node* new_node_pt=0;
7620 if (face0.is_boundary_face())
7622 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7626 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7628 tet_face_node_pt[face0]=new_node_pt;
7629 Node_pt.push_back(new_node_pt);
7630 Vector<double> s(3);
7631 Vector<double> s_tet(3);
7632 Vector<double> x_tet(3);
7633 el_pt->local_coordinate_of_node(j,s);
7634 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7635 tet_el_pt->interpolated_x(s_tet,x_tet);
7636 new_node_pt->x(0)=x_tet[0];
7637 new_node_pt->x(1)=x_tet[1];
7638 new_node_pt->x(2)=x_tet[2];
7643 el_pt->node_pt(j)=old_node_pt;
7652 old_node_pt=tet_face_node_pt[face1];
7655 Node* new_node_pt=0;
7656 if (face1.is_boundary_face())
7658 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7662 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7664 tet_face_node_pt[face1]=new_node_pt;
7665 Node_pt.push_back(new_node_pt);
7666 Vector<double> s(3);
7667 Vector<double> s_tet(3);
7668 Vector<double> x_tet(3);
7669 el_pt->local_coordinate_of_node(j,s);
7670 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7671 tet_el_pt->interpolated_x(s_tet,x_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];
7679 el_pt->node_pt(j)=old_node_pt;
7688 old_node_pt=tet_face_node_pt[face2];
7691 Node* new_node_pt=0;
7692 if (face2.is_boundary_face())
7694 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7698 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7700 tet_face_node_pt[face2]=new_node_pt;
7701 Node_pt.push_back(new_node_pt);
7702 Vector<double> s(3);
7703 Vector<double> s_tet(3);
7704 Vector<double> x_tet(3);
7705 el_pt->local_coordinate_of_node(j,s);
7706 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7707 tet_el_pt->interpolated_x(s_tet,x_tet);
7708 new_node_pt->x(0)=x_tet[0];
7709 new_node_pt->x(1)=x_tet[1];
7710 new_node_pt->x(2)=x_tet[2];
7715 el_pt->node_pt(j)=old_node_pt;
7726 el_pt->node_pt(j)=centroid_node_pt;
7736 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7737 Node_pt.push_back(new_node_pt);
7738 Vector<double> s(3);
7739 Vector<double> s_tet(3);
7740 Vector<double> x_tet(3);
7741 el_pt->local_coordinate_of_node(j,s);
7742 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7743 tet_el_pt->interpolated_x(s_tet,x_tet);
7744 new_node_pt->x(0)=x_tet[0];
7745 new_node_pt->x(1)=x_tet[1];
7746 new_node_pt->x(2)=x_tet[2];
7755 Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
7756 old_node_pt=brick_edge_node_pt[edge];
7759 Node* new_node_pt=0;
7760 if (edge.is_boundary_edge())
7762 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7766 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7768 brick_edge_node_pt[edge]=new_node_pt;
7769 Node_pt.push_back(new_node_pt);
7770 Vector<double> s(3);
7771 Vector<double> s_tet(3);
7772 Vector<double> x_tet(3);
7773 el_pt->local_coordinate_of_node(j,s);
7774 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7775 tet_el_pt->interpolated_x(s_tet,x_tet);
7776 new_node_pt->x(0)=x_tet[0];
7777 new_node_pt->x(1)=x_tet[1];
7778 new_node_pt->x(2)=x_tet[2];
7783 el_pt->node_pt(j)=old_node_pt;
7793 Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
7794 old_node_pt=brick_edge_node_pt[edge];
7797 Node* new_node_pt=0;
7798 if (edge.is_boundary_edge())
7800 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7804 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7806 brick_edge_node_pt[edge]=new_node_pt;
7807 Node_pt.push_back(new_node_pt);
7808 Vector<double> s(3);
7809 Vector<double> s_tet(3);
7810 Vector<double> x_tet(3);
7811 el_pt->local_coordinate_of_node(j,s);
7812 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7813 tet_el_pt->interpolated_x(s_tet,x_tet);
7814 new_node_pt->x(0)=x_tet[0];
7815 new_node_pt->x(1)=x_tet[1];
7816 new_node_pt->x(2)=x_tet[2];
7821 el_pt->node_pt(j)=old_node_pt;
7830 Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
7831 old_node_pt=brick_edge_node_pt[edge];
7834 Node* new_node_pt=0;
7835 if (edge.is_boundary_edge())
7837 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7841 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7843 brick_edge_node_pt[edge]=new_node_pt;
7844 Node_pt.push_back(new_node_pt);
7845 Vector<double> s(3);
7846 Vector<double> s_tet(3);
7847 Vector<double> x_tet(3);
7848 el_pt->local_coordinate_of_node(j,s);
7849 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7850 tet_el_pt->interpolated_x(s_tet,x_tet);
7851 new_node_pt->x(0)=x_tet[0];
7852 new_node_pt->x(1)=x_tet[1];
7853 new_node_pt->x(2)=x_tet[2];
7858 el_pt->node_pt(j)=old_node_pt;
7867 Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
7868 old_node_pt=brick_edge_node_pt[edge];
7871 Node* new_node_pt=0;
7872 if (edge.is_boundary_edge())
7874 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7878 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7880 brick_edge_node_pt[edge]=new_node_pt;
7881 Node_pt.push_back(new_node_pt);
7882 Vector<double> s(3);
7883 Vector<double> s_tet(3);
7884 Vector<double> x_tet(3);
7885 el_pt->local_coordinate_of_node(j,s);
7886 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7887 tet_el_pt->interpolated_x(s_tet,x_tet);
7888 new_node_pt->x(0)=x_tet[0];
7889 new_node_pt->x(1)=x_tet[1];
7890 new_node_pt->x(2)=x_tet[2];
7895 el_pt->node_pt(j)=old_node_pt;
7904 Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
7905 old_node_pt=brick_edge_node_pt[edge];
7908 Node* new_node_pt=0;
7909 if (edge.is_boundary_edge())
7911 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7915 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7917 brick_edge_node_pt[edge]=new_node_pt;
7918 Node_pt.push_back(new_node_pt);
7919 Vector<double> s(3);
7920 Vector<double> s_tet(3);
7921 Vector<double> x_tet(3);
7922 el_pt->local_coordinate_of_node(j,s);
7923 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7924 tet_el_pt->interpolated_x(s_tet,x_tet);
7925 new_node_pt->x(0)=x_tet[0];
7926 new_node_pt->x(1)=x_tet[1];
7927 new_node_pt->x(2)=x_tet[2];
7932 el_pt->node_pt(j)=old_node_pt;
7942 Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
7943 old_node_pt=brick_edge_node_pt[edge];
7946 Node* new_node_pt=0;
7947 if (edge.is_boundary_edge())
7949 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7953 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7955 brick_edge_node_pt[edge]=new_node_pt;
7956 Node_pt.push_back(new_node_pt);
7957 Vector<double> s(3);
7958 Vector<double> s_tet(3);
7959 Vector<double> x_tet(3);
7960 el_pt->local_coordinate_of_node(j,s);
7961 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7962 tet_el_pt->interpolated_x(s_tet,x_tet);
7963 new_node_pt->x(0)=x_tet[0];
7964 new_node_pt->x(1)=x_tet[1];
7965 new_node_pt->x(2)=x_tet[2];
7970 el_pt->node_pt(j)=old_node_pt;
7979 Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
7980 old_node_pt=brick_edge_node_pt[edge];
7983 Node* new_node_pt=0;
7984 if (edge.is_boundary_edge())
7986 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
7990 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
7992 brick_edge_node_pt[edge]=new_node_pt;
7993 Node_pt.push_back(new_node_pt);
7994 Vector<double> s(3);
7995 Vector<double> s_tet(3);
7996 Vector<double> x_tet(3);
7997 el_pt->local_coordinate_of_node(j,s);
7998 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
7999 tet_el_pt->interpolated_x(s_tet,x_tet);
8000 new_node_pt->x(0)=x_tet[0];
8001 new_node_pt->x(1)=x_tet[1];
8002 new_node_pt->x(2)=x_tet[2];
8007 el_pt->node_pt(j)=old_node_pt;
8017 Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
8018 old_node_pt=brick_edge_node_pt[edge];
8021 Node* new_node_pt=0;
8022 if (edge.is_boundary_edge())
8024 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8028 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8030 brick_edge_node_pt[edge]=new_node_pt;
8031 Node_pt.push_back(new_node_pt);
8032 Vector<double> s(3);
8033 Vector<double> s_tet(3);
8034 Vector<double> x_tet(3);
8035 el_pt->local_coordinate_of_node(j,s);
8036 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8037 tet_el_pt->interpolated_x(s_tet,x_tet);
8038 new_node_pt->x(0)=x_tet[0];
8039 new_node_pt->x(1)=x_tet[1];
8040 new_node_pt->x(2)=x_tet[2];
8045 el_pt->node_pt(j)=old_node_pt;
8054 Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
8055 old_node_pt=brick_edge_node_pt[edge];
8058 Node* new_node_pt=0;
8059 if (edge.is_boundary_edge())
8061 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8065 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8067 brick_edge_node_pt[edge]=new_node_pt;
8068 Node_pt.push_back(new_node_pt);
8069 Vector<double> s(3);
8070 Vector<double> s_tet(3);
8071 Vector<double> x_tet(3);
8072 el_pt->local_coordinate_of_node(j,s);
8073 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8074 tet_el_pt->interpolated_x(s_tet,x_tet);
8075 new_node_pt->x(0)=x_tet[0];
8076 new_node_pt->x(1)=x_tet[1];
8077 new_node_pt->x(2)=x_tet[2];
8082 el_pt->node_pt(j)=old_node_pt;
8092 Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
8093 old_node_pt=brick_edge_node_pt[edge];
8096 Node* new_node_pt=0;
8097 if (edge.is_boundary_edge())
8099 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8103 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8105 brick_edge_node_pt[edge]=new_node_pt;
8106 Node_pt.push_back(new_node_pt);
8107 Vector<double> s(3);
8108 Vector<double> s_tet(3);
8109 Vector<double> x_tet(3);
8110 el_pt->local_coordinate_of_node(j,s);
8111 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8112 tet_el_pt->interpolated_x(s_tet,x_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];
8120 el_pt->node_pt(j)=old_node_pt;
8130 Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
8131 old_node_pt=brick_edge_node_pt[edge];
8134 Node* new_node_pt=0;
8135 if (edge.is_boundary_edge())
8137 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8141 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8143 brick_edge_node_pt[edge]=new_node_pt;
8144 Node_pt.push_back(new_node_pt);
8145 Vector<double> s(3);
8146 Vector<double> s_tet(3);
8147 Vector<double> x_tet(3);
8148 el_pt->local_coordinate_of_node(j,s);
8149 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8150 tet_el_pt->interpolated_x(s_tet,x_tet);
8151 new_node_pt->x(0)=x_tet[0];
8152 new_node_pt->x(1)=x_tet[1];
8153 new_node_pt->x(2)=x_tet[2];
8158 el_pt->node_pt(j)=old_node_pt;
8168 Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
8169 old_node_pt=brick_edge_node_pt[edge];
8172 Node* new_node_pt=0;
8173 if (edge.is_boundary_edge())
8175 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8179 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8181 brick_edge_node_pt[edge]=new_node_pt;
8182 Node_pt.push_back(new_node_pt);
8183 Vector<double> s(3);
8184 Vector<double> s_tet(3);
8185 Vector<double> x_tet(3);
8186 el_pt->local_coordinate_of_node(j,s);
8187 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8188 tet_el_pt->interpolated_x(s_tet,x_tet);
8189 new_node_pt->x(0)=x_tet[0];
8190 new_node_pt->x(1)=x_tet[1];
8191 new_node_pt->x(2)=x_tet[2];
8196 el_pt->node_pt(j)=old_node_pt;
8207 TFace face(tet_el_pt->node_pt( central_tet_vertex),
8208 tet_el_pt->node_pt(tet_edge_node[0][0]),
8209 tet_el_pt->node_pt(tet_edge_node[2][0]));
8211 old_node_pt=tet_face_node_pt[face];
8214 Node* new_node_pt=0;
8215 if (face.is_boundary_face())
8217 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8221 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8223 tet_face_node_pt[face]=new_node_pt;
8224 Node_pt.push_back(new_node_pt);
8225 Vector<double> s(3);
8226 Vector<double> s_tet(3);
8227 Vector<double> x_tet(3);
8228 el_pt->local_coordinate_of_node(j,s);
8229 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8230 tet_el_pt->interpolated_x(s_tet,x_tet);
8231 new_node_pt->x(0)=x_tet[0];
8232 new_node_pt->x(1)=x_tet[1];
8233 new_node_pt->x(2)=x_tet[2];
8238 el_pt->node_pt(j)=old_node_pt;
8250 TFace face(tet_el_pt->node_pt(central_tet_vertex),
8251 tet_el_pt->node_pt(tet_edge_node[1][0]),
8252 tet_el_pt->node_pt(tet_edge_node[2][0]));
8253 old_node_pt=tet_face_node_pt[face];
8256 Node* new_node_pt=0;
8257 if (face.is_boundary_face())
8259 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8263 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8265 tet_face_node_pt[face]=new_node_pt;
8266 Node_pt.push_back(new_node_pt);
8267 Vector<double> s(3);
8268 Vector<double> s_tet(3);
8269 Vector<double> x_tet(3);
8270 el_pt->local_coordinate_of_node(j,s);
8271 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8272 tet_el_pt->interpolated_x(s_tet,x_tet);
8273 new_node_pt->x(0)=x_tet[0];
8274 new_node_pt->x(1)=x_tet[1];
8275 new_node_pt->x(2)=x_tet[2];
8280 el_pt->node_pt(j)=old_node_pt;
8292 TFace face(tet_el_pt->node_pt(central_tet_vertex),
8293 tet_el_pt->node_pt(tet_edge_node[0][0]),
8294 tet_el_pt->node_pt(tet_edge_node[1][0]));
8296 old_node_pt=tet_face_node_pt[face];
8299 Node* new_node_pt=0;
8300 if (face.is_boundary_face())
8302 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8306 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8308 tet_face_node_pt[face]=new_node_pt;
8309 Node_pt.push_back(new_node_pt);
8310 Vector<double> s(3);
8311 Vector<double> s_tet(3);
8312 Vector<double> x_tet(3);
8313 el_pt->local_coordinate_of_node(j,s);
8314 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8315 tet_el_pt->interpolated_x(s_tet,x_tet);
8316 new_node_pt->x(0)=x_tet[0];
8317 new_node_pt->x(1)=x_tet[1];
8318 new_node_pt->x(2)=x_tet[2];
8323 el_pt->node_pt(j)=old_node_pt;
8331 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8332 Node_pt.push_back(new_node_pt);
8333 Vector<double> s(3);
8334 Vector<double> s_tet(3);
8335 Vector<double> x_tet(3);
8336 el_pt->local_coordinate_of_node(j,s);
8337 dummy_q_el_pt[2]->interpolated_s_tet(s,s_tet);
8338 top_mid_face_node2_pt=new_node_pt;
8339 tet_el_pt->interpolated_x(s_tet,x_tet);
8340 new_node_pt->x(0)=x_tet[0];
8341 new_node_pt->x(1)=x_tet[1];
8342 new_node_pt->x(2)=x_tet[2];
8351 el_pt->node_pt(j)=back_mid_face_node0_pt;
8360 el_pt->node_pt(j)=right_mid_face_node1_pt;
8370 for (
unsigned j=0;j<8;j++)
8372 Node* nod_pt=dummy_q_el_pt[3]->node_pt(j);
8373 Vector<double> s_tet(3);
8374 Vector<double> x_tet(3);
8378 tet_el_pt->local_coordinate_of_node(2,s_tet);
8379 nod_pt->set_value(0,s_tet[0]);
8380 nod_pt->set_value(1,s_tet[1]);
8381 nod_pt->set_value(2,s_tet[2]);
8382 tet_el_pt->interpolated_x(s_tet,x_tet);
8383 nod_pt->x(0)=x_tet[0];
8384 nod_pt->x(1)=x_tet[1];
8385 nod_pt->x(2)=x_tet[2];
8388 tet_el_pt->local_coordinate_of_node(7,s_tet);
8389 nod_pt->set_value(0,s_tet[0]);
8390 nod_pt->set_value(1,s_tet[1]);
8391 nod_pt->set_value(2,s_tet[2]);
8392 tet_el_pt->interpolated_x(s_tet,x_tet);
8393 nod_pt->x(0)=x_tet[0];
8394 nod_pt->x(1)=x_tet[1];
8395 nod_pt->x(2)=x_tet[2];
8398 tet_el_pt->local_coordinate_of_node(5,s_tet);
8399 nod_pt->set_value(0,s_tet[0]);
8400 nod_pt->set_value(1,s_tet[1]);
8401 nod_pt->set_value(2,s_tet[2]);
8402 tet_el_pt->interpolated_x(s_tet,x_tet);
8403 nod_pt->x(0)=x_tet[0];
8404 nod_pt->x(1)=x_tet[1];
8405 nod_pt->x(2)=x_tet[2];
8413 nod_pt->set_value(0,s_tet[0]);
8414 nod_pt->set_value(1,s_tet[1]);
8415 nod_pt->set_value(2,s_tet[2]);
8416 tet_el_pt->interpolated_x(s_tet,x_tet);
8417 nod_pt->x(0)=x_tet[0];
8418 nod_pt->x(1)=x_tet[1];
8419 nod_pt->x(2)=x_tet[2];
8422 tet_el_pt->local_coordinate_of_node(8,s_tet);
8423 nod_pt->set_value(0,s_tet[0]);
8424 nod_pt->set_value(1,s_tet[1]);
8425 nod_pt->set_value(2,s_tet[2]);
8426 tet_el_pt->interpolated_x(s_tet,x_tet);
8427 nod_pt->x(0)=x_tet[0];
8428 nod_pt->x(1)=x_tet[1];
8429 nod_pt->x(2)=x_tet[2];
8437 nod_pt->set_value(0,s_tet[0]);
8438 nod_pt->set_value(1,s_tet[1]);
8439 nod_pt->set_value(2,s_tet[2]);
8440 tet_el_pt->interpolated_x(s_tet,x_tet);
8441 nod_pt->x(0)=x_tet[0];
8442 nod_pt->x(1)=x_tet[1];
8443 nod_pt->x(2)=x_tet[2];
8451 nod_pt->set_value(0,s_tet[0]);
8452 nod_pt->set_value(1,s_tet[1]);
8453 nod_pt->set_value(2,s_tet[2]);
8454 tet_el_pt->interpolated_x(s_tet,x_tet);
8455 nod_pt->x(0)=x_tet[0];
8456 nod_pt->x(1)=x_tet[1];
8457 nod_pt->x(2)=x_tet[2];
8464 nod_pt->set_value(0,s_tet[0]);
8465 nod_pt->set_value(1,s_tet[1]);
8466 nod_pt->set_value(2,s_tet[2]);
8467 tet_el_pt->interpolated_x(s_tet,x_tet);
8468 nod_pt->x(0)=x_tet[0];
8469 nod_pt->x(1)=x_tet[1];
8470 nod_pt->x(2)=x_tet[2];
8478 FiniteElement* el_pt=
new ELEMENT;
8480 Element_pt.push_back(el_pt);
8482 TFace face0(tet_el_pt->node_pt(1),
8483 tet_el_pt->node_pt(2),
8484 tet_el_pt->node_pt(3));
8486 TFace face1(tet_el_pt->node_pt(0),
8487 tet_el_pt->node_pt(2),
8488 tet_el_pt->node_pt(3));
8490 TFace face2(tet_el_pt->node_pt(0),
8491 tet_el_pt->node_pt(1),
8492 tet_el_pt->node_pt(2));
8495 Vector<Vector<unsigned> > tet_edge_node(3);
8496 tet_edge_node[0].resize(2);
8497 tet_edge_node[0][0]=7;
8498 tet_edge_node[0][1]=1;
8499 tet_edge_node[1].resize(2);
8500 tet_edge_node[1][0]=5;
8501 tet_edge_node[1][1]=0;
8502 tet_edge_node[2].resize(2);
8503 tet_edge_node[2][0]=8;
8504 tet_edge_node[2][1]=3;
8507 unsigned central_tet_vertex=2;
8509 Node* tet_node_pt=0;
8510 Node* old_node_pt=0;
8517 tet_node_pt=tet_el_pt->node_pt(central_tet_vertex);
8518 old_node_pt=tet_node_node_pt[tet_node_pt];
8521 Node* new_node_pt=0;
8522 if (tet_node_pt->is_on_boundary())
8524 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8528 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8530 tet_node_node_pt[tet_node_pt]=new_node_pt;
8531 Node_pt.push_back(new_node_pt);
8532 Vector<double> s(3);
8533 Vector<double> s_tet(3);
8534 Vector<double> x_tet(3);
8535 el_pt->local_coordinate_of_node(j,s);
8536 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8537 tet_el_pt->interpolated_x(s_tet,x_tet);
8538 new_node_pt->x(0)=x_tet[0];
8539 new_node_pt->x(1)=x_tet[1];
8540 new_node_pt->x(2)=x_tet[2];
8545 el_pt->node_pt(j)=old_node_pt;
8555 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[0][0]);
8556 old_node_pt=tet_node_node_pt[tet_node_pt];
8559 Node* new_node_pt=0;
8560 if (tet_node_pt->is_on_boundary())
8562 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8566 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8568 tet_node_node_pt[tet_node_pt]=new_node_pt;
8569 Node_pt.push_back(new_node_pt);
8570 Vector<double> s(3);
8571 Vector<double> s_tet(3);
8572 Vector<double> x_tet(3);
8573 el_pt->local_coordinate_of_node(j,s);
8574 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8575 tet_el_pt->interpolated_x(s_tet,x_tet);
8576 new_node_pt->x(0)=x_tet[0];
8577 new_node_pt->x(1)=x_tet[1];
8578 new_node_pt->x(2)=x_tet[2];
8583 el_pt->node_pt(j)=old_node_pt;
8593 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[1][0]);
8594 old_node_pt=tet_node_node_pt[tet_node_pt];
8597 Node* new_node_pt=0;
8598 if (tet_node_pt->is_on_boundary())
8600 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8604 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8606 tet_node_node_pt[tet_node_pt]=new_node_pt;
8607 Node_pt.push_back(new_node_pt);
8608 Vector<double> s(3);
8609 Vector<double> s_tet(3);
8610 Vector<double> x_tet(3);
8611 el_pt->local_coordinate_of_node(j,s);
8612 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8613 tet_el_pt->interpolated_x(s_tet,x_tet);
8614 new_node_pt->x(0)=x_tet[0];
8615 new_node_pt->x(1)=x_tet[1];
8616 new_node_pt->x(2)=x_tet[2];
8621 el_pt->node_pt(j)=old_node_pt;
8631 tet_node_pt=tet_el_pt->node_pt(tet_edge_node[2][0]);
8632 old_node_pt=tet_node_node_pt[tet_node_pt];
8635 Node* new_node_pt=0;
8636 if (tet_node_pt->is_on_boundary())
8638 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8642 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8644 tet_node_node_pt[tet_node_pt]=new_node_pt;
8645 Node_pt.push_back(new_node_pt);
8646 Vector<double> s(3);
8647 Vector<double> s_tet(3);
8648 Vector<double> x_tet(3);
8649 el_pt->local_coordinate_of_node(j,s);
8650 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8651 tet_el_pt->interpolated_x(s_tet,x_tet);
8652 new_node_pt->x(0)=x_tet[0];
8653 new_node_pt->x(1)=x_tet[1];
8654 new_node_pt->x(2)=x_tet[2];
8659 el_pt->node_pt(j)=old_node_pt;
8670 old_node_pt=tet_face_node_pt[face0];
8673 Node* new_node_pt=0;
8674 if (face0.is_boundary_face())
8676 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8680 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8682 tet_face_node_pt[face0]=new_node_pt;
8683 Node_pt.push_back(new_node_pt);
8684 Vector<double> s(3);
8685 Vector<double> s_tet(3);
8686 Vector<double> x_tet(3);
8687 el_pt->local_coordinate_of_node(j,s);
8688 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8689 tet_el_pt->interpolated_x(s_tet,x_tet);
8690 new_node_pt->x(0)=x_tet[0];
8691 new_node_pt->x(1)=x_tet[1];
8692 new_node_pt->x(2)=x_tet[2];
8697 el_pt->node_pt(j)=old_node_pt;
8706 old_node_pt=tet_face_node_pt[face1];
8709 Node* new_node_pt=0;
8710 if (face1.is_boundary_face())
8712 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8716 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8718 tet_face_node_pt[face1]=new_node_pt;
8719 Node_pt.push_back(new_node_pt);
8720 Vector<double> s(3);
8721 Vector<double> s_tet(3);
8722 Vector<double> x_tet(3);
8723 el_pt->local_coordinate_of_node(j,s);
8724 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8725 tet_el_pt->interpolated_x(s_tet,x_tet);
8726 new_node_pt->x(0)=x_tet[0];
8727 new_node_pt->x(1)=x_tet[1];
8728 new_node_pt->x(2)=x_tet[2];
8733 el_pt->node_pt(j)=old_node_pt;
8742 old_node_pt=tet_face_node_pt[face2];
8745 Node* new_node_pt=0;
8746 if (face2.is_boundary_face())
8748 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8752 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8754 tet_face_node_pt[face2]=new_node_pt;
8755 Node_pt.push_back(new_node_pt);
8756 Vector<double> s(3);
8757 Vector<double> s_tet(3);
8758 Vector<double> x_tet(3);
8759 el_pt->local_coordinate_of_node(j,s);
8760 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8761 tet_el_pt->interpolated_x(s_tet,x_tet);
8762 new_node_pt->x(0)=x_tet[0];
8763 new_node_pt->x(1)=x_tet[1];
8764 new_node_pt->x(2)=x_tet[2];
8769 el_pt->node_pt(j)=old_node_pt;
8780 el_pt->node_pt(j)=centroid_node_pt;
8790 Node* new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8791 Node_pt.push_back(new_node_pt);
8792 Vector<double> s(3);
8793 Vector<double> s_tet(3);
8794 Vector<double> x_tet(3);
8795 el_pt->local_coordinate_of_node(j,s);
8796 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8797 tet_el_pt->interpolated_x(s_tet,x_tet);
8798 new_node_pt->x(0)=x_tet[0];
8799 new_node_pt->x(1)=x_tet[1];
8800 new_node_pt->x(2)=x_tet[2];
8809 Edge edge(el_pt->node_pt(0),el_pt->node_pt(2));
8810 old_node_pt=brick_edge_node_pt[edge];
8813 Node* new_node_pt=0;
8814 if (edge.is_boundary_edge())
8816 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8820 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8822 brick_edge_node_pt[edge]=new_node_pt;
8823 Node_pt.push_back(new_node_pt);
8824 Vector<double> s(3);
8825 Vector<double> s_tet(3);
8826 Vector<double> x_tet(3);
8827 el_pt->local_coordinate_of_node(j,s);
8828 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8829 tet_el_pt->interpolated_x(s_tet,x_tet);
8830 new_node_pt->x(0)=x_tet[0];
8831 new_node_pt->x(1)=x_tet[1];
8832 new_node_pt->x(2)=x_tet[2];
8837 el_pt->node_pt(j)=old_node_pt;
8847 Edge edge(el_pt->node_pt(0),el_pt->node_pt(6));
8848 old_node_pt=brick_edge_node_pt[edge];
8851 Node* new_node_pt=0;
8852 if (edge.is_boundary_edge())
8854 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8858 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8860 brick_edge_node_pt[edge]=new_node_pt;
8861 Node_pt.push_back(new_node_pt);
8862 Vector<double> s(3);
8863 Vector<double> s_tet(3);
8864 Vector<double> x_tet(3);
8865 el_pt->local_coordinate_of_node(j,s);
8866 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8867 tet_el_pt->interpolated_x(s_tet,x_tet);
8868 new_node_pt->x(0)=x_tet[0];
8869 new_node_pt->x(1)=x_tet[1];
8870 new_node_pt->x(2)=x_tet[2];
8875 el_pt->node_pt(j)=old_node_pt;
8884 Edge edge(el_pt->node_pt(2),el_pt->node_pt(8));
8885 old_node_pt=brick_edge_node_pt[edge];
8888 Node* new_node_pt=0;
8889 if (edge.is_boundary_edge())
8891 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8895 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8897 brick_edge_node_pt[edge]=new_node_pt;
8898 Node_pt.push_back(new_node_pt);
8899 Vector<double> s(3);
8900 Vector<double> s_tet(3);
8901 Vector<double> x_tet(3);
8902 el_pt->local_coordinate_of_node(j,s);
8903 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8904 tet_el_pt->interpolated_x(s_tet,x_tet);
8905 new_node_pt->x(0)=x_tet[0];
8906 new_node_pt->x(1)=x_tet[1];
8907 new_node_pt->x(2)=x_tet[2];
8912 el_pt->node_pt(j)=old_node_pt;
8921 Edge edge(el_pt->node_pt(6),el_pt->node_pt(8));
8922 old_node_pt=brick_edge_node_pt[edge];
8925 Node* new_node_pt=0;
8926 if (edge.is_boundary_edge())
8928 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8932 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8934 brick_edge_node_pt[edge]=new_node_pt;
8935 Node_pt.push_back(new_node_pt);
8936 Vector<double> s(3);
8937 Vector<double> s_tet(3);
8938 Vector<double> x_tet(3);
8939 el_pt->local_coordinate_of_node(j,s);
8940 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8941 tet_el_pt->interpolated_x(s_tet,x_tet);
8942 new_node_pt->x(0)=x_tet[0];
8943 new_node_pt->x(1)=x_tet[1];
8944 new_node_pt->x(2)=x_tet[2];
8949 el_pt->node_pt(j)=old_node_pt;
8958 Edge edge(el_pt->node_pt(18),el_pt->node_pt(20));
8959 old_node_pt=brick_edge_node_pt[edge];
8962 Node* new_node_pt=0;
8963 if (edge.is_boundary_edge())
8965 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
8969 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
8971 brick_edge_node_pt[edge]=new_node_pt;
8972 Node_pt.push_back(new_node_pt);
8973 Vector<double> s(3);
8974 Vector<double> s_tet(3);
8975 Vector<double> x_tet(3);
8976 el_pt->local_coordinate_of_node(j,s);
8977 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
8978 tet_el_pt->interpolated_x(s_tet,x_tet);
8979 new_node_pt->x(0)=x_tet[0];
8980 new_node_pt->x(1)=x_tet[1];
8981 new_node_pt->x(2)=x_tet[2];
8986 el_pt->node_pt(j)=old_node_pt;
8996 Edge edge(el_pt->node_pt(18),el_pt->node_pt(24));
8997 old_node_pt=brick_edge_node_pt[edge];
9000 Node* new_node_pt=0;
9001 if (edge.is_boundary_edge())
9003 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9007 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9009 brick_edge_node_pt[edge]=new_node_pt;
9010 Node_pt.push_back(new_node_pt);
9011 Vector<double> s(3);
9012 Vector<double> s_tet(3);
9013 Vector<double> x_tet(3);
9014 el_pt->local_coordinate_of_node(j,s);
9015 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9016 tet_el_pt->interpolated_x(s_tet,x_tet);
9017 new_node_pt->x(0)=x_tet[0];
9018 new_node_pt->x(1)=x_tet[1];
9019 new_node_pt->x(2)=x_tet[2];
9024 el_pt->node_pt(j)=old_node_pt;
9033 Edge edge(el_pt->node_pt(20),el_pt->node_pt(26));
9034 old_node_pt=brick_edge_node_pt[edge];
9037 Node* new_node_pt=0;
9038 if (edge.is_boundary_edge())
9040 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9044 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9046 brick_edge_node_pt[edge]=new_node_pt;
9047 Node_pt.push_back(new_node_pt);
9048 Vector<double> s(3);
9049 Vector<double> s_tet(3);
9050 Vector<double> x_tet(3);
9051 el_pt->local_coordinate_of_node(j,s);
9052 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9053 tet_el_pt->interpolated_x(s_tet,x_tet);
9054 new_node_pt->x(0)=x_tet[0];
9055 new_node_pt->x(1)=x_tet[1];
9056 new_node_pt->x(2)=x_tet[2];
9061 el_pt->node_pt(j)=old_node_pt;
9071 Edge edge(el_pt->node_pt(24),el_pt->node_pt(26));
9072 old_node_pt=brick_edge_node_pt[edge];
9075 Node* new_node_pt=0;
9076 if (edge.is_boundary_edge())
9078 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9082 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9084 brick_edge_node_pt[edge]=new_node_pt;
9085 Node_pt.push_back(new_node_pt);
9086 Vector<double> s(3);
9087 Vector<double> s_tet(3);
9088 Vector<double> x_tet(3);
9089 el_pt->local_coordinate_of_node(j,s);
9090 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9091 tet_el_pt->interpolated_x(s_tet,x_tet);
9092 new_node_pt->x(0)=x_tet[0];
9093 new_node_pt->x(1)=x_tet[1];
9094 new_node_pt->x(2)=x_tet[2];
9099 el_pt->node_pt(j)=old_node_pt;
9108 Edge edge(el_pt->node_pt(0),el_pt->node_pt(18));
9109 old_node_pt=brick_edge_node_pt[edge];
9112 Node* new_node_pt=0;
9113 if (edge.is_boundary_edge())
9115 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9119 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9121 brick_edge_node_pt[edge]=new_node_pt;
9122 Node_pt.push_back(new_node_pt);
9123 Vector<double> s(3);
9124 Vector<double> s_tet(3);
9125 Vector<double> x_tet(3);
9126 el_pt->local_coordinate_of_node(j,s);
9127 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9128 tet_el_pt->interpolated_x(s_tet,x_tet);
9129 new_node_pt->x(0)=x_tet[0];
9130 new_node_pt->x(1)=x_tet[1];
9131 new_node_pt->x(2)=x_tet[2];
9136 el_pt->node_pt(j)=old_node_pt;
9146 Edge edge(el_pt->node_pt(2),el_pt->node_pt(20));
9147 old_node_pt=brick_edge_node_pt[edge];
9150 Node* new_node_pt=0;
9151 if (edge.is_boundary_edge())
9153 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9157 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9159 brick_edge_node_pt[edge]=new_node_pt;
9160 Node_pt.push_back(new_node_pt);
9161 Vector<double> s(3);
9162 Vector<double> s_tet(3);
9163 Vector<double> x_tet(3);
9164 el_pt->local_coordinate_of_node(j,s);
9165 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9166 tet_el_pt->interpolated_x(s_tet,x_tet);
9167 new_node_pt->x(0)=x_tet[0];
9168 new_node_pt->x(1)=x_tet[1];
9169 new_node_pt->x(2)=x_tet[2];
9174 el_pt->node_pt(j)=old_node_pt;
9184 Edge edge(el_pt->node_pt(6),el_pt->node_pt(24));
9185 old_node_pt=brick_edge_node_pt[edge];
9188 Node* new_node_pt=0;
9189 if (edge.is_boundary_edge())
9191 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9195 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9197 brick_edge_node_pt[edge]=new_node_pt;
9198 Node_pt.push_back(new_node_pt);
9199 Vector<double> s(3);
9200 Vector<double> s_tet(3);
9201 Vector<double> x_tet(3);
9202 el_pt->local_coordinate_of_node(j,s);
9203 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9204 tet_el_pt->interpolated_x(s_tet,x_tet);
9205 new_node_pt->x(0)=x_tet[0];
9206 new_node_pt->x(1)=x_tet[1];
9207 new_node_pt->x(2)=x_tet[2];
9212 el_pt->node_pt(j)=old_node_pt;
9222 Edge edge(el_pt->node_pt(8),el_pt->node_pt(26));
9223 old_node_pt=brick_edge_node_pt[edge];
9226 Node* new_node_pt=0;
9227 if (edge.is_boundary_edge())
9229 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9233 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9235 brick_edge_node_pt[edge]=new_node_pt;
9236 Node_pt.push_back(new_node_pt);
9237 Vector<double> s(3);
9238 Vector<double> s_tet(3);
9239 Vector<double> x_tet(3);
9240 el_pt->local_coordinate_of_node(j,s);
9241 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9242 tet_el_pt->interpolated_x(s_tet,x_tet);
9243 new_node_pt->x(0)=x_tet[0];
9244 new_node_pt->x(1)=x_tet[1];
9245 new_node_pt->x(2)=x_tet[2];
9250 el_pt->node_pt(j)=old_node_pt;
9261 TFace face(tet_el_pt->node_pt( central_tet_vertex),
9262 tet_el_pt->node_pt(tet_edge_node[0][0]),
9263 tet_el_pt->node_pt(tet_edge_node[2][0]));
9265 old_node_pt=tet_face_node_pt[face];
9268 Node* new_node_pt=0;
9269 if (face.is_boundary_face())
9271 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9275 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9277 tet_face_node_pt[face]=new_node_pt;
9278 Node_pt.push_back(new_node_pt);
9279 Vector<double> s(3);
9280 Vector<double> s_tet(3);
9281 Vector<double> x_tet(3);
9282 el_pt->local_coordinate_of_node(j,s);
9283 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9284 tet_el_pt->interpolated_x(s_tet,x_tet);
9285 new_node_pt->x(0)=x_tet[0];
9286 new_node_pt->x(1)=x_tet[1];
9287 new_node_pt->x(2)=x_tet[2];
9292 el_pt->node_pt(j)=old_node_pt;
9304 TFace face(tet_el_pt->node_pt(central_tet_vertex),
9305 tet_el_pt->node_pt(tet_edge_node[1][0]),
9306 tet_el_pt->node_pt(tet_edge_node[2][0]));
9308 old_node_pt=tet_face_node_pt[face];
9311 Node* new_node_pt=0;
9312 if (face.is_boundary_face())
9314 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9318 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9320 tet_face_node_pt[face]=new_node_pt;
9321 Node_pt.push_back(new_node_pt);
9322 Vector<double> s(3);
9323 Vector<double> s_tet(3);
9324 Vector<double> x_tet(3);
9325 el_pt->local_coordinate_of_node(j,s);
9326 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9327 tet_el_pt->interpolated_x(s_tet,x_tet);
9328 new_node_pt->x(0)=x_tet[0];
9329 new_node_pt->x(1)=x_tet[1];
9330 new_node_pt->x(2)=x_tet[2];
9335 el_pt->node_pt(j)=old_node_pt;
9347 TFace face(tet_el_pt->node_pt(central_tet_vertex),
9348 tet_el_pt->node_pt(tet_edge_node[0][0]),
9349 tet_el_pt->node_pt(tet_edge_node[1][0]));
9351 old_node_pt=tet_face_node_pt[face];
9354 Node* new_node_pt=0;
9355 if (face.is_boundary_face())
9357 new_node_pt=el_pt->construct_boundary_node(j,time_stepper_pt);
9361 new_node_pt=el_pt->construct_node(j,time_stepper_pt);
9363 tet_face_node_pt[face]=new_node_pt;
9364 Node_pt.push_back(new_node_pt);
9365 Vector<double> s(3);
9366 Vector<double> s_tet(3);
9367 Vector<double> x_tet(3);
9368 el_pt->local_coordinate_of_node(j,s);
9369 dummy_q_el_pt[3]->interpolated_s_tet(s,s_tet);
9370 tet_el_pt->interpolated_x(s_tet,x_tet);
9371 new_node_pt->x(0)=x_tet[0];
9372 new_node_pt->x(1)=x_tet[1];
9373 new_node_pt->x(2)=x_tet[2];
9378 el_pt->node_pt(j)=old_node_pt;
9388 el_pt->node_pt(j)=top_mid_face_node2_pt;
9397 el_pt->node_pt(j)=top_mid_face_node1_pt;
9406 el_pt->node_pt(j)=top_mid_face_node0_pt;
9418 for (
int face_index=0;face_index<4;face_index++)
9428 face_pt=
new TFace(tet_el_pt->node_pt(1),
9429 tet_el_pt->node_pt(2),
9430 tet_el_pt->node_pt(3));
9435 face_pt=
new TFace(tet_el_pt->node_pt(0),
9436 tet_el_pt->node_pt(2),
9437 tet_el_pt->node_pt(3));
9444 face_pt=
new TFace(tet_el_pt->node_pt(0),
9445 tet_el_pt->node_pt(1),
9446 tet_el_pt->node_pt(3));
9452 face_pt=
new TFace(tet_el_pt->node_pt(0),
9453 tet_el_pt->node_pt(1),
9454 tet_el_pt->node_pt(2));
9460 if (face_pt->is_boundary_face())
9462 std::set<unsigned> bnd;
9463 std::set<unsigned>* bnd_pt=&bnd;
9464 face_pt->get_boundaries_pt(bnd_pt);
9467 if ((*bnd_pt).size()>1)
9469 std::ostringstream error_stream;
9470 error_stream <<
"TFace should only be on one boundary.\n";
9471 throw OomphLibError(
9473 OOMPH_CURRENT_FUNCTION,
9474 OOMPH_EXCEPTION_LOCATION);
9478 if ((*bnd_pt).size()==1)
9481 FaceElement* face_el_pt=0;
9482 if (tet_mesh_is_solid_mesh)
9485 if (
dynamic_cast<SolidTElement<3,3>*
>(tet_el_pt)==0)
9487 std::ostringstream error_stream;
9489 <<
"Tet-element cannot be cast to SolidTElement<3,3>.\n" 9490 <<
"BrickFromTetMesh can only be built from\n" 9491 <<
"mesh containing quadratic tets.\n" 9493 throw OomphLibError(error_stream.str(),
9494 OOMPH_CURRENT_FUNCTION,
9495 OOMPH_EXCEPTION_LOCATION);
9500 new DummyFaceElement<SolidTElement<3,3> >(tet_el_pt,face_index);
9505 if (
dynamic_cast<TElement<3,3>*
>(tet_el_pt)==0)
9507 std::ostringstream error_stream;
9509 <<
"Tet-element cannot be cast to TElement<3,3>.\n" 9510 <<
"BrickFromTetMesh can only be built from\n" 9511 <<
"mesh containing quadratic tets.\n" 9513 throw OomphLibError(error_stream.str(),
9514 OOMPH_CURRENT_FUNCTION,
9515 OOMPH_EXCEPTION_LOCATION);
9521 new DummyFaceElement<TElement<3,3> >(tet_el_pt,face_index);
9527 unsigned b=(*(*bnd_pt).begin());
9528 Boundary_coordinate_exists[b]=
true;
9529 face_el_pt->set_boundary_number_in_bulk_mesh(b);
9533 Vector<Node*> brick_face_node_pt(19);
9538 brick_face_node_pt[0]=brick_el1_pt->node_pt(0);
9539 brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
9540 brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
9542 brick_face_node_pt[3]=brick_el1_pt->node_pt(18);
9543 brick_face_node_pt[4]=brick_el2_pt->node_pt(18);
9544 brick_face_node_pt[5]=brick_el1_pt->node_pt(2);
9546 brick_face_node_pt[6]=brick_el1_pt->node_pt(9);
9547 brick_face_node_pt[7]=brick_el3_pt->node_pt(1);
9548 brick_face_node_pt[8]=brick_el3_pt->node_pt(9);
9550 brick_face_node_pt[9]=brick_el2_pt->node_pt(9);
9551 brick_face_node_pt[10]=brick_el2_pt->node_pt(3);
9552 brick_face_node_pt[11]=brick_el1_pt->node_pt(1);
9554 brick_face_node_pt[12]=brick_el1_pt->node_pt(20);
9556 brick_face_node_pt[13]=brick_el2_pt->node_pt(12);
9557 brick_face_node_pt[14]=brick_el1_pt->node_pt(19);
9559 brick_face_node_pt[15]=brick_el1_pt->node_pt(10);
9560 brick_face_node_pt[16]=brick_el2_pt->node_pt(21);
9562 brick_face_node_pt[17]=brick_el3_pt->node_pt(10);
9563 brick_face_node_pt[18]=brick_el1_pt->node_pt(11);
9567 brick_face_node_pt[0]=brick_el0_pt->node_pt(0);
9568 brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
9569 brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
9571 brick_face_node_pt[3]=brick_el0_pt->node_pt(18);
9572 brick_face_node_pt[4]=brick_el2_pt->node_pt(18);
9573 brick_face_node_pt[5]=brick_el0_pt->node_pt(6);
9575 brick_face_node_pt[6]=brick_el0_pt->node_pt(9);
9576 brick_face_node_pt[7]=brick_el3_pt->node_pt(3);
9577 brick_face_node_pt[8]=brick_el3_pt->node_pt(9);
9579 brick_face_node_pt[9]=brick_el2_pt->node_pt(9);
9580 brick_face_node_pt[10]=brick_el2_pt->node_pt(1);
9581 brick_face_node_pt[11]=brick_el0_pt->node_pt(3);
9583 brick_face_node_pt[12]=brick_el0_pt->node_pt(24);
9585 brick_face_node_pt[13]=brick_el2_pt->node_pt(10);
9586 brick_face_node_pt[14]=brick_el0_pt->node_pt(21);
9588 brick_face_node_pt[15]=brick_el0_pt->node_pt(12);
9589 brick_face_node_pt[16]=brick_el3_pt->node_pt(21);
9591 brick_face_node_pt[17]=brick_el3_pt->node_pt(12);
9592 brick_face_node_pt[18]=brick_el0_pt->node_pt(15);
9596 brick_face_node_pt[0]=brick_el0_pt->node_pt(0);
9597 brick_face_node_pt[1]=brick_el1_pt->node_pt(0);
9598 brick_face_node_pt[2]=brick_el2_pt->node_pt(0);
9600 brick_face_node_pt[3]=brick_el0_pt->node_pt(2);
9601 brick_face_node_pt[4]=brick_el1_pt->node_pt(2);
9602 brick_face_node_pt[5]=brick_el0_pt->node_pt(6);
9604 brick_face_node_pt[6]=brick_el0_pt->node_pt(1);
9605 brick_face_node_pt[7]=brick_el1_pt->node_pt(3);
9606 brick_face_node_pt[8]=brick_el1_pt->node_pt(1);
9608 brick_face_node_pt[9]=brick_el2_pt->node_pt(3);
9609 brick_face_node_pt[10]=brick_el2_pt->node_pt(1);
9610 brick_face_node_pt[11]=brick_el0_pt->node_pt(3);
9612 brick_face_node_pt[12]=brick_el0_pt->node_pt(8);
9614 brick_face_node_pt[13]=brick_el2_pt->node_pt(4);
9615 brick_face_node_pt[14]=brick_el0_pt->node_pt(5);
9617 brick_face_node_pt[15]=brick_el0_pt->node_pt(4);
9618 brick_face_node_pt[16]=brick_el1_pt->node_pt(5);
9620 brick_face_node_pt[17]=brick_el1_pt->node_pt(4);
9621 brick_face_node_pt[18]=brick_el0_pt->node_pt(7);
9625 brick_face_node_pt[0]=brick_el1_pt->node_pt(0);
9626 brick_face_node_pt[1]=brick_el3_pt->node_pt(0);
9627 brick_face_node_pt[2]=brick_el0_pt->node_pt(0);
9629 brick_face_node_pt[3]=brick_el1_pt->node_pt(18);
9630 brick_face_node_pt[4]=brick_el3_pt->node_pt(6);
9631 brick_face_node_pt[5]=brick_el1_pt->node_pt(6);
9633 brick_face_node_pt[6]=brick_el1_pt->node_pt(9);
9634 brick_face_node_pt[7]=brick_el3_pt->node_pt(1);
9635 brick_face_node_pt[8]=brick_el3_pt->node_pt(3);
9637 brick_face_node_pt[9]=brick_el0_pt->node_pt(9);
9638 brick_face_node_pt[10]=brick_el0_pt->node_pt(1);
9639 brick_face_node_pt[11]=brick_el1_pt->node_pt(3);
9641 brick_face_node_pt[12]=brick_el1_pt->node_pt(24);
9643 brick_face_node_pt[13]=brick_el0_pt->node_pt(10);
9644 brick_face_node_pt[14]=brick_el1_pt->node_pt(21);
9646 brick_face_node_pt[15]=brick_el1_pt->node_pt(12);
9647 brick_face_node_pt[16]=brick_el3_pt->node_pt(7);
9649 brick_face_node_pt[17]=brick_el3_pt->node_pt(4);
9650 brick_face_node_pt[18]=brick_el1_pt->node_pt(15);
9657 Vector<unsigned> translate(19);
9660 for (
unsigned i=0;i<19;i++)
9666 for (
unsigned j=0;j<19;j++)
9669 Node* brick_node_pt=brick_face_node_pt[translate[j]];
9672 Vector<double> s=s_face[j];
9673 Vector<double> zeta(2);
9674 Vector<double> x(3);
9675 face_el_pt->interpolated_zeta(s,zeta);
9676 face_el_pt->interpolated_x(s,x);
9680 double dist=sqrt(pow(brick_node_pt->x(0)-x[0],2)+
9681 pow(brick_node_pt->x(1)-x[1],2)+
9682 pow(brick_node_pt->x(2)-x[2],2));
9683 if (dist>BrickFromTetMeshHelper::Face_position_tolerance)
9685 std::ofstream brick0;
9686 std::ofstream brick1;
9687 std::ofstream brick2;
9688 std::ofstream brick3;
9689 brick0.open(
"full_brick0.dat");
9690 brick1.open(
"full_brick1.dat");
9691 brick2.open(
"full_brick2.dat");
9692 brick3.open(
"full_brick3.dat");
9693 for (
unsigned j=0;j<27;j++)
9695 brick0 << brick_el0_pt->node_pt(j)->x(0) <<
" " 9696 << brick_el0_pt->node_pt(j)->x(1) <<
" " 9697 << brick_el0_pt->node_pt(j)->x(2) <<
"\n";
9699 brick1 << brick_el1_pt->node_pt(j)->x(0) <<
" " 9700 << brick_el1_pt->node_pt(j)->x(1) <<
" " 9701 << brick_el1_pt->node_pt(j)->x(2) <<
"\n";
9703 brick2 << brick_el2_pt->node_pt(j)->x(0) <<
" " 9704 << brick_el2_pt->node_pt(j)->x(1) <<
" " 9705 << brick_el2_pt->node_pt(j)->x(2) <<
"\n";
9707 brick3 << brick_el3_pt->node_pt(j)->x(0) <<
" " 9708 << brick_el3_pt->node_pt(j)->x(1) <<
" " 9709 << brick_el3_pt->node_pt(j)->x(2) <<
"\n";
9716 std::ofstream full_face;
9717 full_face.open(
"full_face.dat");
9718 for (
unsigned j=0;j<6;j++)
9720 full_face << face_el_pt->node_pt(j)->x(0) <<
" " 9721 << face_el_pt->node_pt(j)->x(1) <<
" " 9722 << face_el_pt->node_pt(j)->x(2) <<
"\n";
9727 int normal_sign=face_el_pt->normal_sign();
9729 std::ostringstream error_stream;
9731 <<
"During assignment of boundary cordinates, the distance\n" 9732 <<
"between brick node and reference point in \n" 9733 <<
"triangular FaceElement is " << dist <<
" which \n" 9734 <<
"is bigger than the tolerance defined in \n" 9735 <<
"BrickFromTetMeshHelper::Face_position_tolerance=" 9736 << BrickFromTetMeshHelper::Face_position_tolerance <<
".\n" 9737 <<
"If this is tolerable, increase the tolerance \n" 9738 <<
"(it's defined in a namespace and therefore publically\n" 9739 <<
"accessible). If not, the Face may be inverted in which \n" 9740 <<
"case you should re-implement the translation scheme,\n" 9741 <<
"following the pattern used in the ThinLayerBrickOnTetMesh." 9742 <<
"\nThe required code fragements are already here but \n" 9743 <<
"the translation is the unit map.\n" 9744 <<
"To aid the diagnostics, the files full_brick[0-3].dat\n" 9745 <<
"contain the coordinates of the 27 nodes in the four\n" 9746 <<
"bricks associated with the current tet and full_face.dat\n" 9747 <<
"contains the coordinates of the 6 nodes in the FaceElement" 9748 <<
"\nfrom which the boundary coordinates are extracted.\n" 9749 <<
"FYI: The normal_sign of the face is: " << normal_sign
9751 throw OomphLibError(
9753 OOMPH_CURRENT_FUNCTION,
9754 OOMPH_EXCEPTION_LOCATION);
9759 add_boundary_node(b,brick_node_pt);
9760 brick_node_pt->set_coordinates_on_boundary(b,zeta);
9767 Boundary_element_pt[b].push_back(brick_el1_pt);
9768 Face_index_at_boundary[b].push_back(-2);
9769 Boundary_element_pt[b].push_back(brick_el2_pt);
9770 Face_index_at_boundary[b].push_back(-1);
9771 Boundary_element_pt[b].push_back(brick_el3_pt);
9772 Face_index_at_boundary[b].push_back(-2);
9776 Boundary_element_pt[b].push_back(brick_el0_pt);
9777 Face_index_at_boundary[b].push_back(-1);
9778 Boundary_element_pt[b].push_back(brick_el2_pt);
9779 Face_index_at_boundary[b].push_back(-2);
9780 Boundary_element_pt[b].push_back(brick_el3_pt);
9781 Face_index_at_boundary[b].push_back(-1);
9785 Boundary_element_pt[b].push_back(brick_el0_pt);
9786 Face_index_at_boundary[b].push_back(-3);
9787 Boundary_element_pt[b].push_back(brick_el1_pt);
9788 Face_index_at_boundary[b].push_back(-3);
9789 Boundary_element_pt[b].push_back(brick_el2_pt);
9790 Face_index_at_boundary[b].push_back(-3);
9794 Boundary_element_pt[b].push_back(brick_el0_pt);
9795 Face_index_at_boundary[b].push_back(-2);
9796 Boundary_element_pt[b].push_back(brick_el1_pt);
9797 Face_index_at_boundary[b].push_back(-1);
9798 Boundary_element_pt[b].push_back(brick_el3_pt);
9799 Face_index_at_boundary[b].push_back(-3);
9812 Lookup_for_elements_next_boundary_is_setup=
true;
9828 for (
unsigned e=0;e<4;e++)
9830 for (
unsigned j=0;j<8;j++)
9832 delete dummy_q_el_pt[e]->node_pt(j);
9834 delete dummy_q_el_pt[e];
Unstructured tet mesh based on output from Tetgen: http://wias-berlin.de/software/tetgen/.
void build_mesh(XdaTetMesh< TElement< 3, 3 > > *tet_mesh_pt, TimeStepper *time_stepper_pt)
Build fct: Pass pointer to existing tet mesh.
Tet mesh made of quadratic (ten node) tets built from xda input file.