26 #ifndef OOMPH_THIN_LAYER_BRICK_ON_TET_MESH_TEMPLATE_CC
27 #define OOMPH_THIN_LAYER_BRICK_ON_TET_MESH_TEMPLATE_CC
30 #include "../solid/solid_elements.h"
52 template<
class ELEMENT>
55 const Vector<unsigned>& boundary_ids,
56 ThicknessFctPt thickness_fct_pt,
57 const unsigned& nlayer,
58 const Vector<Vector<unsigned>>& in_out_boundary_id,
59 TimeStepper* time_stepper_pt)
60 : Thickness_fct_pt(thickness_fct_pt)
63 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
66 bool tet_mesh_is_solid_mesh =
false;
67 if (
dynamic_cast<SolidFiniteElement*
>(tet_mesh_pt->element_pt(0)) != 0)
69 tet_mesh_is_solid_mesh =
true;
76 Vector<Vector<double>> s_face(19);
77 for (
unsigned i = 0; i < 19; i++)
154 s_face[i][0] = 1.0 / 3.0;
155 s_face[i][1] = 1.0 / 3.0;
162 s_face[i][0] = 5.0 / 24.0;
163 s_face[i][1] = 5.0 / 24.0;
167 s_face[i][0] = 5.0 / 12.0;
168 s_face[i][1] = 5.0 / 12.0;
174 s_face[i][1] = 5.0 / 24.0;
175 s_face[i][0] = 7.0 / 12.0;
179 s_face[i][1] = 5.0 / 12.0;
180 s_face[i][0] = 1.0 / 6.0;
187 s_face[i][0] = 5.0 / 24.0;
188 s_face[i][1] = 7.0 / 12.0;
192 s_face[i][0] = 5.0 / 12.0;
193 s_face[i][1] = 1.0 / 6.0;
200 MapMatrixMixed<int, unsigned, unsigned> translate;
203 for (
unsigned i = 0; i < 19; i++)
205 translate(-1, i) = i;
208 translate(-1, 6) = 11;
209 translate(-1, 11) = 6;
210 translate(-1, 3) = 5;
211 translate(-1, 5) = 3;
212 translate(-1, 18) = 14;
213 translate(-1, 14) = 18;
214 translate(-1, 7) = 10;
215 translate(-1, 10) = 7;
216 translate(-1, 13) = 17;
217 translate(-1, 17) = 13;
218 translate(-1, 1) = 2;
219 translate(-1, 2) = 1;
220 translate(-1, 9) = 8;
221 translate(-1, 8) = 9;
225 std::map<Node*, Node*> solid_node_pt;
228 std::map<Edge, Node*> quarter_edge_node;
232 std::map<Node*, Vector<Vector<double>>> normals;
236 std::map<Node*, Vector<Node*>> connected_node_pt;
239 Element_pt.reserve(3 * boundary_ids.size() * nlayer);
243 std::set<unsigned> all_bnd;
246 unsigned nb = boundary_ids.size();
247 for (
unsigned ib = 0; ib < nb; ib++)
250 unsigned b = boundary_ids[ib];
254 unsigned nnod = tet_mesh_pt->nboundary_node(b);
255 for (
unsigned j = 0; j < nnod; j++)
257 Node* nod_pt = tet_mesh_pt->boundary_node_pt(b, j);
260 std::set<unsigned>* bnd_pt;
261 nod_pt->get_boundaries_pt(bnd_pt);
264 for (std::set<unsigned>::iterator it = (*bnd_pt).begin();
265 it != (*bnd_pt).end();
274 unsigned highest_fluid_bound_id =
275 *std::max_element(all_bnd.begin(), all_bnd.end());
278 std::vector<bool> is_on_fsi_boundary(highest_fluid_bound_id + 1,
false);
279 for (
unsigned ib = 0; ib < nb; ib++)
281 is_on_fsi_boundary[boundary_ids[ib]] =
true;
287 Vector<std::vector<bool>> is_on_in_out_boundary(n);
288 Vector<std::set<unsigned>> in_out_boundary_id_set(n);
289 for (
unsigned j = 0; j < n; j++)
291 is_on_in_out_boundary[j].resize(highest_fluid_bound_id + 1,
false);
293 for (
unsigned ib = 0; ib < nb; ib++)
305 unsigned maxb = highest_fluid_bound_id + 2;
311 Boundary_element_pt.resize(maxb);
312 Face_index_at_boundary.resize(maxb);
316 nb = boundary_ids.size();
317 for (
unsigned ib = 0; ib < nb; ib++)
320 unsigned b = boundary_ids[ib];
324 Boundary_coordinate_exists[b] =
true;
330 unsigned nel = tet_mesh_pt->nboundary_element(b);
331 for (
unsigned e = 0; e < nel; e++)
334 FiniteElement* bulk_elem_pt = tet_mesh_pt->boundary_element_pt(b, e);
337 int face_index = tet_mesh_pt->face_index_at_boundary(b, e);
340 FaceElement* face_el_pt = 0;
341 if (tet_mesh_is_solid_mesh)
344 if (
dynamic_cast<SolidTElement<3, 3>*
>(bulk_elem_pt) == 0)
346 std::ostringstream error_stream;
348 <<
"Tet-element cannot be cast to SolidTElement<3,3>.\n"
349 <<
"ThinBrickOnTetMesh can only be erected on mesh containing\n"
350 <<
"quadratic tets." << std::endl;
351 throw OomphLibError(error_stream.str(),
352 OOMPH_CURRENT_FUNCTION,
353 OOMPH_EXCEPTION_LOCATION);
357 new DummyFaceElement<SolidTElement<3, 3>>(bulk_elem_pt, face_index);
362 if (
dynamic_cast<TElement<3, 3>*
>(bulk_elem_pt) == 0)
364 std::ostringstream error_stream;
366 <<
"Tet-element cannot be cast to TElement<3,3>.\n"
367 <<
"ThinBrickOnTetMesh can only be erected on mesh containing\n"
368 <<
"quadratic tets." << std::endl;
369 throw OomphLibError(error_stream.str(),
370 OOMPH_CURRENT_FUNCTION,
371 OOMPH_EXCEPTION_LOCATION);
375 new DummyFaceElement<TElement<3, 3>>(bulk_elem_pt, face_index);
381 face_el_pt->set_boundary_number_in_bulk_mesh(b);
384 Vector<Vector<FiniteElement*>> new_el_pt(3);
393 for (
unsigned j = 0; j < 3; j++)
396 new_el_pt[j].resize(nlayer);
397 for (
unsigned ilayer = 0; ilayer < nlayer; ilayer++)
399 new_el_pt[j][ilayer] =
new ELEMENT;
400 Element_pt.push_back(new_el_pt[j][ilayer]);
403 Boundary_element_pt[b].push_back(new_el_pt[j][0]);
404 Face_index_at_boundary[b].push_back(-3);
408 unsigned j_local = 0;
411 normal_sign = face_el_pt->normal_sign();
415 Vector<double> s = s_face[translate(normal_sign, j)];
416 Vector<double> zeta(2);
418 Vector<double> unit_normal(3);
419 face_el_pt->interpolated_zeta(s, zeta);
420 face_el_pt->interpolated_x(s, x);
421 face_el_pt->outer_unit_normal(s, unit_normal);
425 Node* fluid_node_pt = face_el_pt->node_pt(translate(normal_sign, j));
428 Node* existing_node_pt = solid_node_pt[fluid_node_pt];
429 if (existing_node_pt == 0)
432 Node* new_node_pt = new_el_pt[j][0]->construct_boundary_node(
433 j_local, time_stepper_pt);
434 Node_pt.push_back(new_node_pt);
437 solid_node_pt[fluid_node_pt] = new_node_pt;
440 new_node_pt->x(0) = x[0];
441 new_node_pt->x(1) = x[1];
442 new_node_pt->x(2) = x[2];
445 bool only_on_fsi =
true;
446 std::set<unsigned>* bnd_pt;
447 fluid_node_pt->get_boundaries_pt(bnd_pt);
448 for (std::set<unsigned>::iterator it = (*bnd_pt).begin();
449 it != (*bnd_pt).end();
452 if (!is_on_fsi_boundary[(*it)]) only_on_fsi =
false;
453 add_boundary_node((*it), new_node_pt);
455 new_node_pt->set_coordinates_on_boundary(b, zeta);
456 normals[new_node_pt].push_back(unit_normal);
465 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
466 connected_node_pt[new_node_pt].push_back(new_nod_pt);
467 Node_pt.push_back(new_nod_pt);
472 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
473 j_local + 18, time_stepper_pt);
474 connected_node_pt[new_node_pt].push_back(new_nod_pt);
475 Node_pt.push_back(new_nod_pt);
479 Node* new_nod_pt = new_el_pt[j][0]->construct_node(
480 j_local + 18, time_stepper_pt);
481 connected_node_pt[new_node_pt].push_back(new_nod_pt);
482 Node_pt.push_back(new_nod_pt);
486 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
489 new_el_pt[j][ilayer]->node_pt(j_local) =
490 connected_node_pt[new_node_pt][2 * ilayer - 1];
493 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
494 j_local + 9, time_stepper_pt);
495 connected_node_pt[new_node_pt].push_back(new_nod_pt);
496 Node_pt.push_back(new_nod_pt);
499 if (ilayer != (nlayer - 1))
501 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
502 j_local + 18, time_stepper_pt);
503 connected_node_pt[new_node_pt].push_back(new_nod_pt);
504 Node_pt.push_back(new_nod_pt);
509 new_el_pt[j][ilayer]->construct_boundary_node(
510 j_local + 18, time_stepper_pt);
511 connected_node_pt[new_node_pt].push_back(new_nod_pt);
512 Node_pt.push_back(new_nod_pt);
519 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
520 j_local + 9, time_stepper_pt);
521 connected_node_pt[new_node_pt].push_back(new_nod_pt);
522 Node_pt.push_back(new_nod_pt);
524 new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
525 j_local + 18, time_stepper_pt);
526 connected_node_pt[new_node_pt].push_back(new_nod_pt);
527 Node_pt.push_back(new_nod_pt);
530 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
533 new_el_pt[j][ilayer]->node_pt(j_local) =
534 connected_node_pt[new_node_pt][2 * ilayer - 1];
538 new_el_pt[j][ilayer]->construct_boundary_node(
539 j_local + 9, time_stepper_pt);
540 connected_node_pt[new_node_pt].push_back(new_nod_pt);
541 Node_pt.push_back(new_nod_pt);
542 new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
543 j_local + 18, time_stepper_pt);
544 connected_node_pt[new_node_pt].push_back(new_nod_pt);
545 Node_pt.push_back(new_nod_pt);
553 existing_node_pt->set_coordinates_on_boundary(b, zeta);
554 normals[existing_node_pt].push_back(unit_normal);
557 new_el_pt[j][0]->node_pt(j_local) = existing_node_pt;
558 new_el_pt[j][0]->node_pt(j_local + 9) =
559 connected_node_pt[existing_node_pt][0];
560 new_el_pt[j][0]->node_pt(j_local + 18) =
561 connected_node_pt[existing_node_pt][1];
564 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
566 new_el_pt[j][ilayer]->node_pt(j_local) =
567 connected_node_pt[existing_node_pt][2 * ilayer - 1];
568 new_el_pt[j][ilayer]->node_pt(j_local + 9) =
569 connected_node_pt[existing_node_pt][2 * ilayer];
570 new_el_pt[j][ilayer]->node_pt(j_local + 18) =
571 connected_node_pt[existing_node_pt][2 * ilayer + 1];
582 s = s_face[translate(normal_sign, j + 3)];
583 face_el_pt->interpolated_zeta(s, zeta);
584 face_el_pt->interpolated_x(s, x);
585 face_el_pt->outer_unit_normal(s, unit_normal);
588 fluid_node_pt = face_el_pt->node_pt(translate(normal_sign, j + 3));
591 existing_node_pt = solid_node_pt[fluid_node_pt];
592 if (existing_node_pt == 0)
595 Node* new_node_pt = new_el_pt[j][0]->construct_boundary_node(
596 j_local, time_stepper_pt);
597 Node_pt.push_back(new_node_pt);
600 solid_node_pt[fluid_node_pt] = new_node_pt;
603 new_node_pt->x(0) = x[0];
604 new_node_pt->x(1) = x[1];
605 new_node_pt->x(2) = x[2];
608 bool only_on_fsi =
true;
609 std::set<unsigned>* bnd_pt;
610 fluid_node_pt->get_boundaries_pt(bnd_pt);
611 for (std::set<unsigned>::iterator it = (*bnd_pt).begin();
612 it != (*bnd_pt).end();
615 if (!is_on_fsi_boundary[(*it)]) only_on_fsi =
false;
616 add_boundary_node((*it), new_node_pt);
618 new_node_pt->set_coordinates_on_boundary(b, zeta);
619 normals[new_node_pt].push_back(unit_normal);
627 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
628 connected_node_pt[new_node_pt].push_back(new_nod_pt);
629 Node_pt.push_back(new_nod_pt);
634 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
635 j_local + 18, time_stepper_pt);
636 connected_node_pt[new_node_pt].push_back(new_nod_pt);
637 Node_pt.push_back(new_nod_pt);
641 Node* new_nod_pt = new_el_pt[j][0]->construct_node(
642 j_local + 18, time_stepper_pt);
643 connected_node_pt[new_node_pt].push_back(new_nod_pt);
644 Node_pt.push_back(new_nod_pt);
648 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
651 new_el_pt[j][ilayer]->node_pt(j_local) =
652 connected_node_pt[new_node_pt][2 * ilayer - 1];
655 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
656 j_local + 9, time_stepper_pt);
657 connected_node_pt[new_node_pt].push_back(new_nod_pt);
658 Node_pt.push_back(new_nod_pt);
661 if (ilayer != (nlayer - 1))
663 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
664 j_local + 18, time_stepper_pt);
665 connected_node_pt[new_node_pt].push_back(new_nod_pt);
666 Node_pt.push_back(new_nod_pt);
671 new_el_pt[j][ilayer]->construct_boundary_node(
672 j_local + 18, time_stepper_pt);
673 connected_node_pt[new_node_pt].push_back(new_nod_pt);
674 Node_pt.push_back(new_nod_pt);
681 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
682 j_local + 9, time_stepper_pt);
683 connected_node_pt[new_node_pt].push_back(new_nod_pt);
684 Node_pt.push_back(new_nod_pt);
686 new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
687 j_local + 18, time_stepper_pt);
688 connected_node_pt[new_node_pt].push_back(new_nod_pt);
689 Node_pt.push_back(new_nod_pt);
692 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
695 new_el_pt[j][ilayer]->node_pt(j_local) =
696 connected_node_pt[new_node_pt][2 * ilayer - 1];
700 new_el_pt[j][ilayer]->construct_boundary_node(
701 j_local + 9, time_stepper_pt);
702 connected_node_pt[new_node_pt].push_back(new_nod_pt);
703 Node_pt.push_back(new_nod_pt);
704 new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
705 j_local + 18, time_stepper_pt);
706 connected_node_pt[new_node_pt].push_back(new_nod_pt);
707 Node_pt.push_back(new_nod_pt);
715 existing_node_pt->set_coordinates_on_boundary(b, zeta);
716 normals[existing_node_pt].push_back(unit_normal);
719 new_el_pt[j][0]->node_pt(j_local) = existing_node_pt;
720 new_el_pt[j][0]->node_pt(j_local + 9) =
721 connected_node_pt[existing_node_pt][0];
722 new_el_pt[j][0]->node_pt(j_local + 18) =
723 connected_node_pt[existing_node_pt][1];
726 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
728 new_el_pt[j][ilayer]->node_pt(j_local) =
729 connected_node_pt[existing_node_pt][2 * ilayer - 1];
730 new_el_pt[j][ilayer]->node_pt(j_local + 9) =
731 connected_node_pt[existing_node_pt][2 * ilayer];
732 new_el_pt[j][ilayer]->node_pt(j_local + 18) =
733 connected_node_pt[existing_node_pt][2 * ilayer + 1];
744 s = s_face[translate(normal_sign, 6 + 2 * j)];
745 face_el_pt->interpolated_zeta(s, zeta);
746 face_el_pt->interpolated_x(s, x);
747 face_el_pt->outer_unit_normal(s, unit_normal);
750 Edge edge(face_el_pt->node_pt(translate(normal_sign, j)),
751 face_el_pt->node_pt(translate(normal_sign, j + 3)));
754 existing_node_pt = quarter_edge_node[edge];
755 if (existing_node_pt == 0)
758 Node* new_node_pt = new_el_pt[j][0]->construct_boundary_node(
759 j_local, time_stepper_pt);
760 Node_pt.push_back(new_node_pt);
763 quarter_edge_node[edge] = new_node_pt;
766 new_node_pt->x(0) = x[0];
767 new_node_pt->x(1) = x[1];
768 new_node_pt->x(2) = x[2];
771 std::set<unsigned>* bnd1_pt;
772 edge.node1_pt()->get_boundaries_pt(bnd1_pt);
773 std::set<unsigned>* bnd2_pt;
774 edge.node2_pt()->get_boundaries_pt(bnd2_pt);
775 std::set<unsigned> bnd;
776 set_intersection((*bnd1_pt).begin(),
780 inserter(bnd, bnd.begin()));
781 bool only_on_fsi =
true;
782 for (std::set<unsigned>::iterator it = bnd.begin(); it != bnd.end();
785 if (!is_on_fsi_boundary[(*it)]) only_on_fsi =
false;
786 add_boundary_node((*it), new_node_pt);
788 new_node_pt->set_coordinates_on_boundary(b, zeta);
789 normals[new_node_pt].push_back(unit_normal);
798 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
799 connected_node_pt[new_node_pt].push_back(new_nod_pt);
800 Node_pt.push_back(new_nod_pt);
805 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
806 j_local + 18, time_stepper_pt);
807 connected_node_pt[new_node_pt].push_back(new_nod_pt);
808 Node_pt.push_back(new_nod_pt);
812 Node* new_nod_pt = new_el_pt[j][0]->construct_node(
813 j_local + 18, time_stepper_pt);
814 connected_node_pt[new_node_pt].push_back(new_nod_pt);
815 Node_pt.push_back(new_nod_pt);
819 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
822 new_el_pt[j][ilayer]->node_pt(j_local) =
823 connected_node_pt[new_node_pt][2 * ilayer - 1];
826 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
827 j_local + 9, time_stepper_pt);
828 connected_node_pt[new_node_pt].push_back(new_nod_pt);
829 Node_pt.push_back(new_nod_pt);
832 if (ilayer != (nlayer - 1))
834 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
835 j_local + 18, time_stepper_pt);
836 connected_node_pt[new_node_pt].push_back(new_nod_pt);
837 Node_pt.push_back(new_nod_pt);
842 new_el_pt[j][ilayer]->construct_boundary_node(
843 j_local + 18, time_stepper_pt);
844 connected_node_pt[new_node_pt].push_back(new_nod_pt);
845 Node_pt.push_back(new_nod_pt);
852 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
853 j_local + 9, time_stepper_pt);
854 connected_node_pt[new_node_pt].push_back(new_nod_pt);
855 Node_pt.push_back(new_nod_pt);
857 new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
858 j_local + 18, time_stepper_pt);
859 connected_node_pt[new_node_pt].push_back(new_nod_pt);
860 Node_pt.push_back(new_nod_pt);
863 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
866 new_el_pt[j][ilayer]->node_pt(j_local) =
867 connected_node_pt[new_node_pt][2 * ilayer - 1];
871 new_el_pt[j][ilayer]->construct_boundary_node(
872 j_local + 9, time_stepper_pt);
873 connected_node_pt[new_node_pt].push_back(new_nod_pt);
874 Node_pt.push_back(new_nod_pt);
875 new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
876 j_local + 18, time_stepper_pt);
877 connected_node_pt[new_node_pt].push_back(new_nod_pt);
878 Node_pt.push_back(new_nod_pt);
886 existing_node_pt->set_coordinates_on_boundary(b, zeta);
887 normals[existing_node_pt].push_back(unit_normal);
890 new_el_pt[j][0]->node_pt(j_local) = existing_node_pt;
891 new_el_pt[j][0]->node_pt(j_local + 9) =
892 connected_node_pt[existing_node_pt][0];
893 new_el_pt[j][0]->node_pt(j_local + 18) =
894 connected_node_pt[existing_node_pt][1];
898 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
900 new_el_pt[j][ilayer]->node_pt(j_local) =
901 connected_node_pt[existing_node_pt][2 * ilayer - 1];
902 new_el_pt[j][ilayer]->node_pt(j_local + 9) =
903 connected_node_pt[existing_node_pt][2 * ilayer];
904 new_el_pt[j][ilayer]->node_pt(j_local + 18) =
905 connected_node_pt[existing_node_pt][2 * ilayer + 1];
915 unsigned other_node = 0;
932 Edge edge2(face_el_pt->node_pt(translate(normal_sign, j)),
933 face_el_pt->node_pt(translate(normal_sign, other_node)));
936 s = s_face[translate(normal_sign, jj)];
937 face_el_pt->interpolated_zeta(s, zeta);
938 face_el_pt->interpolated_x(s, x);
939 face_el_pt->outer_unit_normal(s, unit_normal);
942 existing_node_pt = quarter_edge_node[edge2];
943 if (existing_node_pt == 0)
946 Node* new_node_pt = new_el_pt[j][0]->construct_boundary_node(
947 j_local, time_stepper_pt);
948 Node_pt.push_back(new_node_pt);
951 quarter_edge_node[edge2] = new_node_pt;
954 new_node_pt->x(0) = x[0];
955 new_node_pt->x(1) = x[1];
956 new_node_pt->x(2) = x[2];
959 std::set<unsigned>* bnd1_pt;
960 edge2.node1_pt()->get_boundaries_pt(bnd1_pt);
961 std::set<unsigned>* bnd2_pt;
962 edge2.node2_pt()->get_boundaries_pt(bnd2_pt);
963 std::set<unsigned> bnd;
964 set_intersection((*bnd1_pt).begin(),
968 inserter(bnd, bnd.begin()));
969 bool only_on_fsi =
true;
970 for (std::set<unsigned>::iterator it = bnd.begin(); it != bnd.end();
973 if (!is_on_fsi_boundary[(*it)]) only_on_fsi =
false;
974 add_boundary_node((*it), new_node_pt);
976 new_node_pt->set_coordinates_on_boundary(b, zeta);
977 normals[new_node_pt].push_back(unit_normal);
985 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
986 connected_node_pt[new_node_pt].push_back(new_nod_pt);
987 Node_pt.push_back(new_nod_pt);
992 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
993 j_local + 18, time_stepper_pt);
994 connected_node_pt[new_node_pt].push_back(new_nod_pt);
995 Node_pt.push_back(new_nod_pt);
999 Node* new_nod_pt = new_el_pt[j][0]->construct_node(
1000 j_local + 18, time_stepper_pt);
1001 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1002 Node_pt.push_back(new_nod_pt);
1006 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1009 new_el_pt[j][ilayer]->node_pt(j_local) =
1010 connected_node_pt[new_node_pt][2 * ilayer - 1];
1013 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1014 j_local + 9, time_stepper_pt);
1015 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1016 Node_pt.push_back(new_nod_pt);
1018 if (ilayer != (nlayer - 1))
1020 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1021 j_local + 18, time_stepper_pt);
1022 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1023 Node_pt.push_back(new_nod_pt);
1028 new_el_pt[j][ilayer]->construct_boundary_node(
1029 j_local + 18, time_stepper_pt);
1030 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1031 Node_pt.push_back(new_nod_pt);
1038 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
1039 j_local + 9, time_stepper_pt);
1040 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1041 Node_pt.push_back(new_nod_pt);
1043 new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
1044 j_local + 18, time_stepper_pt);
1045 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1046 Node_pt.push_back(new_nod_pt);
1049 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1052 new_el_pt[j][ilayer]->node_pt(j_local) =
1053 connected_node_pt[new_node_pt][2 * ilayer - 1];
1057 new_el_pt[j][ilayer]->construct_boundary_node(
1058 j_local + 9, time_stepper_pt);
1059 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1060 Node_pt.push_back(new_nod_pt);
1061 new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
1062 j_local + 18, time_stepper_pt);
1063 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1064 Node_pt.push_back(new_nod_pt);
1072 existing_node_pt->set_coordinates_on_boundary(b, zeta);
1073 normals[existing_node_pt].push_back(unit_normal);
1076 new_el_pt[j][0]->node_pt(j_local) = existing_node_pt;
1077 new_el_pt[j][0]->node_pt(j_local + 9) =
1078 connected_node_pt[existing_node_pt][0];
1079 new_el_pt[j][0]->node_pt(j_local + 18) =
1080 connected_node_pt[existing_node_pt][1];
1083 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1085 new_el_pt[j][ilayer]->node_pt(j_local) =
1086 connected_node_pt[existing_node_pt][2 * ilayer - 1];
1087 new_el_pt[j][ilayer]->node_pt(j_local + 9) =
1088 connected_node_pt[existing_node_pt][2 * ilayer];
1089 new_el_pt[j][ilayer]->node_pt(j_local + 18) =
1090 connected_node_pt[existing_node_pt][2 * ilayer + 1];
1101 new_el_pt[j][0]->construct_boundary_node(j_local, time_stepper_pt);
1102 Node_pt.push_back(new_node_pt);
1119 s = s_face[translate(normal_sign, jj)];
1120 face_el_pt->interpolated_zeta(s, zeta);
1121 face_el_pt->interpolated_x(s, x);
1122 face_el_pt->outer_unit_normal(s, unit_normal);
1125 new_node_pt->x(0) = x[0];
1126 new_node_pt->x(1) = x[1];
1127 new_node_pt->x(2) = x[2];
1130 add_boundary_node(b, new_node_pt);
1131 new_node_pt->set_coordinates_on_boundary(b, zeta);
1132 normals[new_node_pt].push_back(unit_normal);
1136 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
1137 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1138 Node_pt.push_back(new_nod_pt);
1143 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
1144 j_local + 18, time_stepper_pt);
1145 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1146 Node_pt.push_back(new_nod_pt);
1151 new_el_pt[j][0]->construct_node(j_local + 18, time_stepper_pt);
1152 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1153 Node_pt.push_back(new_nod_pt);
1157 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1160 new_el_pt[j][ilayer]->node_pt(j_local) =
1161 connected_node_pt[new_node_pt][2 * ilayer - 1];
1164 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1165 j_local + 9, time_stepper_pt);
1166 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1167 Node_pt.push_back(new_nod_pt);
1169 if (ilayer != (nlayer - 1))
1171 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1172 j_local + 18, time_stepper_pt);
1173 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1174 Node_pt.push_back(new_nod_pt);
1178 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
1179 j_local + 18, time_stepper_pt);
1180 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1181 Node_pt.push_back(new_nod_pt);
1195 new_el_pt[j][0]->construct_boundary_node(j_local, time_stepper_pt);
1196 Node_pt.push_back(new_node_pt);
1214 s = s_face[translate(normal_sign, jj)];
1215 face_el_pt->interpolated_zeta(s, zeta);
1216 face_el_pt->interpolated_x(s, x);
1217 face_el_pt->outer_unit_normal(s, unit_normal);
1220 new_node_pt->x(0) = x[0];
1221 new_node_pt->x(1) = x[1];
1222 new_node_pt->x(2) = x[2];
1225 add_boundary_node(b, new_node_pt);
1226 new_node_pt->set_coordinates_on_boundary(b, zeta);
1227 normals[new_node_pt].push_back(unit_normal);
1231 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
1232 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1233 Node_pt.push_back(new_nod_pt);
1238 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
1239 j_local + 18, time_stepper_pt);
1240 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1241 Node_pt.push_back(new_nod_pt);
1246 new_el_pt[j][0]->construct_node(j_local + 18, time_stepper_pt);
1247 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1248 Node_pt.push_back(new_nod_pt);
1252 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1255 new_el_pt[j][ilayer]->node_pt(j_local) =
1256 connected_node_pt[new_node_pt][2 * ilayer - 1];
1259 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1260 j_local + 9, time_stepper_pt);
1261 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1262 Node_pt.push_back(new_nod_pt);
1265 if (ilayer != (nlayer - 1))
1267 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1268 j_local + 18, time_stepper_pt);
1269 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1270 Node_pt.push_back(new_nod_pt);
1274 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
1275 j_local + 18, time_stepper_pt);
1276 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1277 Node_pt.push_back(new_nod_pt);
1287 unsigned j_local = 8;
1291 new_el_pt[2][0]->construct_boundary_node(j_local, time_stepper_pt);
1292 Node_pt.push_back(new_node_pt);
1296 Vector<double> s = s_face[12];
1297 Vector<double> zeta(2);
1298 Vector<double> x(3);
1299 Vector<double> unit_normal(3);
1300 face_el_pt->interpolated_zeta(s, zeta);
1301 face_el_pt->interpolated_x(s, x);
1302 face_el_pt->outer_unit_normal(s, unit_normal);
1305 new_node_pt->x(0) = x[0];
1306 new_node_pt->x(1) = x[1];
1307 new_node_pt->x(2) = x[2];
1310 add_boundary_node(b, new_node_pt);
1311 new_node_pt->set_coordinates_on_boundary(b, zeta);
1312 normals[new_node_pt].push_back(unit_normal);
1316 new_el_pt[2][0]->construct_node(j_local + 9, time_stepper_pt);
1317 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1318 Node_pt.push_back(new_nod_pt);
1323 Node* new_nod_pt = new_el_pt[2][0]->construct_boundary_node(
1324 j_local + 18, time_stepper_pt);
1325 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1326 Node_pt.push_back(new_nod_pt);
1331 new_el_pt[2][0]->construct_node(j_local + 18, time_stepper_pt);
1332 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1333 Node_pt.push_back(new_nod_pt);
1337 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1340 new_el_pt[2][ilayer]->node_pt(j_local) =
1341 connected_node_pt[new_node_pt][2 * ilayer - 1];
1345 new_el_pt[2][ilayer]->construct_node(j_local + 9, time_stepper_pt);
1346 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1347 Node_pt.push_back(new_nod_pt);
1350 if (ilayer != (nlayer - 1))
1352 Node* new_nod_pt = new_el_pt[2][ilayer]->construct_node(
1353 j_local + 18, time_stepper_pt);
1354 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1355 Node_pt.push_back(new_nod_pt);
1359 Node* new_nod_pt = new_el_pt[2][ilayer]->construct_boundary_node(
1360 j_local + 18, time_stepper_pt);
1361 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1362 Node_pt.push_back(new_nod_pt);
1367 new_el_pt[1][0]->node_pt(j_local) = new_node_pt;
1368 new_el_pt[0][0]->node_pt(j_local) = new_node_pt;
1370 new_el_pt[1][0]->node_pt(j_local + 9) =
1371 connected_node_pt[new_node_pt][0];
1372 new_el_pt[0][0]->node_pt(j_local + 9) =
1373 connected_node_pt[new_node_pt][0];
1375 new_el_pt[1][0]->node_pt(j_local + 18) =
1376 connected_node_pt[new_node_pt][1];
1377 new_el_pt[0][0]->node_pt(j_local + 18) =
1378 connected_node_pt[new_node_pt][1];
1381 for (
unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1383 new_el_pt[1][ilayer]->node_pt(j_local) =
1384 connected_node_pt[new_node_pt][2 * ilayer - 1];
1385 new_el_pt[0][ilayer]->node_pt(j_local) =
1386 connected_node_pt[new_node_pt][2 * ilayer - 1];
1388 new_el_pt[1][ilayer]->node_pt(j_local + 9) =
1389 connected_node_pt[new_node_pt][2 * ilayer];
1390 new_el_pt[0][ilayer]->node_pt(j_local + 9) =
1391 connected_node_pt[new_node_pt][2 * ilayer];
1393 new_el_pt[1][ilayer]->node_pt(j_local + 18) =
1394 connected_node_pt[new_node_pt][2 * ilayer + 1];
1395 new_el_pt[0][ilayer]->node_pt(j_local + 18) =
1396 connected_node_pt[new_node_pt][2 * ilayer + 1];
1404 for (
unsigned ilayer = 0; ilayer < nlayer; ilayer++)
1406 for (
unsigned j = 0; j < 3; j++)
1408 unsigned offset = 9 * j;
1409 new_el_pt[2][ilayer]->node_pt(6 + offset) =
1410 new_el_pt[1][ilayer]->node_pt(2 + offset);
1412 new_el_pt[1][ilayer]->node_pt(6 + offset) =
1413 new_el_pt[0][ilayer]->node_pt(2 + offset);
1415 new_el_pt[0][ilayer]->node_pt(6 + offset) =
1416 new_el_pt[2][ilayer]->node_pt(2 + offset);
1418 new_el_pt[2][ilayer]->node_pt(7 + offset) =
1419 new_el_pt[1][ilayer]->node_pt(5 + offset);
1421 new_el_pt[1][ilayer]->node_pt(7 + offset) =
1422 new_el_pt[0][ilayer]->node_pt(5 + offset);
1424 new_el_pt[0][ilayer]->node_pt(7 + offset) =
1425 new_el_pt[2][ilayer]->node_pt(5 + offset);
1434 unsigned nb_in_out = is_on_in_out_boundary.size();
1440 for (
unsigned j_stack = 0; j_stack < 3; j_stack++)
1443 FiniteElement* el_pt = new_el_pt[j_stack][0];
1446 for (
unsigned j = 0; j < 9; j++)
1449 Node* nod_pt = el_pt->node_pt(j);
1450 Vector<Node*> layer_node_pt = connected_node_pt[nod_pt];
1451 unsigned n = layer_node_pt.size();
1454 std::set<unsigned>* bnd_pt;
1455 nod_pt->get_boundaries_pt(bnd_pt);
1458 for (std::set<unsigned>::iterator it = (*bnd_pt).begin();
1459 it != (*bnd_pt).end();
1463 if (!is_on_fsi_boundary[(*it)])
1466 unsigned ilayer = 0;
1467 for (
unsigned k = 0; k < n; k++)
1470 add_boundary_node((*it), layer_node_pt[k]);
1476 if (j == 1) face_index = -2;
1477 if (j == 3) face_index = -1;
1478 if (j == 5) face_index = 1;
1479 if (j == 7) face_index = 2;
1481 if (face_index != 0)
1487 Boundary_element_pt[(*it)].push_back(
1488 new_el_pt[j_stack][ilayer]);
1489 Face_index_at_boundary[(*it)].push_back(face_index);
1499 for (
unsigned jj = 0; jj < nb_in_out; jj++)
1501 if (is_on_in_out_boundary[jj][(*it)])
1503 in_out_boundary_id_set[jj].insert((*it));
1518 new_el_pt[j_stack][nlayer - 1]);
1532 for (
unsigned jj = 0; jj < n; jj++)
1534 for (std::set<unsigned>::iterator it = in_out_boundary_id_set[jj].begin();
1535 it != in_out_boundary_id_set[jj].end();
1545 unsigned nel = Element_pt.size();
1546 for (
unsigned e = 0; e < nel; e++)
1548 FiniteElement* el_pt = finite_element_pt(e);
1549 for (
unsigned j = 0; j < 27; j++)
1551 if (el_pt->node_pt(j) == 0)
1554 std::ostringstream error_stream;
1555 error_stream <<
"Null node in element " << e <<
" node " << j
1557 throw OomphLibError(error_stream.str(),
1558 OOMPH_CURRENT_FUNCTION,
1559 OOMPH_EXCEPTION_LOCATION);
1567 std::ofstream outfile;
1568 bool doc_normals =
false;
1569 if (doc_normals) outfile.open(
"normals.dat");
1570 for (std::map<Node*, Vector<Vector<double>>>::iterator it = normals.begin();
1571 it != normals.end();
1574 Vector<double> unit_normal(3, 0.0);
1575 unsigned nnormals = ((*it).second).size();
1576 for (
unsigned j = 0; j < nnormals; j++)
1578 for (
unsigned i = 0; i < 3; i++)
1580 unit_normal[i] += ((*it).second)[j][i];
1584 for (
unsigned i = 0; i < 3; i++)
1586 norm += unit_normal[i] * unit_normal[i];
1588 for (
unsigned i = 0; i < 3; i++)
1590 unit_normal[i] /= sqrt(norm);
1593 Node* base_node_pt = (*it).first;
1594 Vector<double> base_pos(3);
1595 base_node_pt->position(base_pos);
1598 Vector<Node*> layer_node_pt = connected_node_pt[base_node_pt];
1599 unsigned n = layer_node_pt.size();
1600 for (
unsigned j = 0; j < n; j++)
1602 for (
unsigned i = 0; i < 3; i++)
1604 layer_node_pt[j]->x(i) =
1605 base_pos[i] + h_thick * double(j + 1) / double(n) * unit_normal[i];
1610 outfile << ((*it).first)->x(0) <<
" " << ((*it).first)->x(1) <<
" "
1611 << ((*it).first)->x(2) <<
" " << unit_normal[0] <<
" "
1612 << unit_normal[1] <<
" " << unit_normal[2] <<
"\n";
1615 if (doc_normals) outfile.close();
1619 Lookup_for_elements_next_boundary_is_setup =
true;
Vector< unsigned > in_out_boundary_id(const unsigned &boundary_id)
Access function to the vector containing the ids of the oomph-lib mesh boundaries that make up the sp...
ThicknessFctPt Thickness_fct_pt
Function pointer to function that specifies the wall thickness as a fct of the coordinates of the inn...
unsigned Outer_boundary_id
Boundary ID of the "outer" surface – the non-wetted tube surface at a distance h_thick from the FSI s...
Vector< Vector< unsigned > > In_out_boundary_id
Vector of vectors containing the ids of the oomph-lib mesh boundaries that make up the specified in/o...
ThinLayerBrickOnTetMesh(Mesh *tet_mesh_pt, const Vector< unsigned > &boundary_ids, ThicknessFctPt thickness_fct_pt, const unsigned &nlayer, const Vector< Vector< unsigned >> &in_out_boundary_id, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Specify (quadratic) tet mesh, boundary IDs of boundary on which the current mesh is to b...
Vector< unsigned > FSI_boundary_id
Vector of oomph-lib boundary ids that make up boundary on which the mesh was erected (typically the F...
////////////////////////////////////////////////////////////////////// //////////////////////////////...