Build fct: Pass pointer to existing tet mesh.
Build fct: Pass pointer to existing tet mesh and timestepper Specialisation for TetgenMesh<TElement<3,3> >
4900 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3, 3);
4903 bool tet_mesh_is_solid_mesh =
false;
4904 if (
dynamic_cast<SolidFiniteElement*
>(tet_mesh_pt->element_pt(0)) != 0)
4906 tet_mesh_is_solid_mesh =
true;
4913 Vector<Vector<double>> s_face(19);
4914 for (
unsigned i = 0;
i < 19;
i++)
4916 s_face[
i].resize(2);
4958 s_face[
i][0] = 0.75;
4959 s_face[
i][1] = 0.25;
4963 s_face[
i][0] = 0.25;
4964 s_face[
i][1] = 0.75;
4969 s_face[
i][1] = 0.75;
4974 s_face[
i][1] = 0.25;
4978 s_face[
i][0] = 0.25;
4983 s_face[
i][0] = 0.75;
4990 s_face[
i][0] = 1.0 / 3.0;
4991 s_face[
i][1] = 1.0 / 3.0;
4998 s_face[
i][0] = 5.0 / 24.0;
4999 s_face[
i][1] = 5.0 / 24.0;
5003 s_face[
i][0] = 5.0 / 12.0;
5004 s_face[
i][1] = 5.0 / 12.0;
5010 s_face[
i][1] = 5.0 / 24.0;
5011 s_face[
i][0] = 7.0 / 12.0;
5015 s_face[
i][1] = 5.0 / 12.0;
5016 s_face[
i][0] = 1.0 / 6.0;
5023 s_face[
i][0] = 5.0 / 24.0;
5024 s_face[
i][1] = 7.0 / 12.0;
5028 s_face[
i][0] = 5.0 / 12.0;
5029 s_face[
i][1] = 1.0 / 6.0;
5035 unsigned nb = tet_mesh_pt->nboundary();
5045 std::map<Node*, Node*> tet_node_node_pt;
5048 std::map<Edge, Node*> brick_edge_node_pt;
5051 std::map<TFace, Node*> tet_face_node_pt;
5055 Vector<DummyBrickElement*> dummy_q_el_pt(4);
5056 for (
unsigned e = 0;
e < 4;
e++)
5058 dummy_q_el_pt[
e] =
new DummyBrickElement;
5059 for (
unsigned j = 0;
j < 8;
j++)
5061 dummy_q_el_pt[
e]->construct_node(
j);
5066 unsigned n_el_tet = tet_mesh_pt->nelement();
5067 for (
unsigned e_tet = 0; e_tet < n_el_tet; e_tet++)
5070 TElement<3, 3>* tet_el_pt =
5071 dynamic_cast<TElement<3, 3>*
>(tet_mesh_pt->element_pt(e_tet));
5076 std::ostringstream error_stream;
5078 <<
"BrickFromTetMesh can only built from tet mesh containing\n"
5079 <<
"ten-noded tets.\n";
5080 throw OomphLibError(
5086 Node* centroid_node_pt = 0;
5089 Node* top_mid_face_node0_pt = 0;
5090 Node* right_mid_face_node0_pt = 0;
5091 Node* back_mid_face_node0_pt = 0;
5093 Node* top_mid_face_node1_pt = 0;
5094 Node* right_mid_face_node1_pt = 0;
5096 Node* top_mid_face_node2_pt = 0;
5099 FiniteElement* brick_el0_pt = 0;
5100 FiniteElement* brick_el1_pt = 0;
5101 FiniteElement* brick_el2_pt = 0;
5102 FiniteElement* brick_el3_pt = 0;
5109 for (
unsigned j = 0;
j < 8;
j++)
5111 Node* nod_pt = dummy_q_el_pt[0]->node_pt(
j);
5112 Vector<double> s_tet(3);
5113 Vector<double> x_tet(3);
5117 tet_el_pt->local_coordinate_of_node(0, s_tet);
5118 nod_pt->set_value(0, s_tet[0]);
5119 nod_pt->set_value(1, s_tet[1]);
5120 nod_pt->set_value(2, s_tet[2]);
5121 tet_el_pt->interpolated_x(s_tet, x_tet);
5122 nod_pt->x(0) = x_tet[0];
5123 nod_pt->x(1) = x_tet[1];
5124 nod_pt->x(2) = x_tet[2];
5127 tet_el_pt->local_coordinate_of_node(4, s_tet);
5128 nod_pt->set_value(0, s_tet[0]);
5129 nod_pt->set_value(1, s_tet[1]);
5130 nod_pt->set_value(2, s_tet[2]);
5131 tet_el_pt->interpolated_x(s_tet, x_tet);
5132 nod_pt->x(0) = x_tet[0];
5133 nod_pt->x(1) = x_tet[1];
5134 nod_pt->x(2) = x_tet[2];
5137 tet_el_pt->local_coordinate_of_node(6, s_tet);
5138 nod_pt->set_value(0, s_tet[0]);
5139 nod_pt->set_value(1, s_tet[1]);
5140 nod_pt->set_value(2, s_tet[2]);
5141 tet_el_pt->interpolated_x(s_tet, x_tet);
5142 nod_pt->x(0) = x_tet[0];
5143 nod_pt->x(1) = x_tet[1];
5144 nod_pt->x(2) = x_tet[2];
5149 s_tet[0] = 1.0 / 3.0;
5150 s_tet[1] = 1.0 / 3.0;
5152 nod_pt->set_value(0, s_tet[0]);
5153 nod_pt->set_value(1, s_tet[1]);
5154 nod_pt->set_value(2, s_tet[2]);
5155 tet_el_pt->interpolated_x(s_tet, x_tet);
5156 nod_pt->x(0) = x_tet[0];
5157 nod_pt->x(1) = x_tet[1];
5158 nod_pt->x(2) = x_tet[2];
5161 tet_el_pt->local_coordinate_of_node(5, s_tet);
5162 nod_pt->set_value(0, s_tet[0]);
5163 nod_pt->set_value(1, s_tet[1]);
5164 nod_pt->set_value(2, s_tet[2]);
5165 tet_el_pt->interpolated_x(s_tet, x_tet);
5166 nod_pt->x(0) = x_tet[0];
5167 nod_pt->x(1) = x_tet[1];
5168 nod_pt->x(2) = x_tet[2];
5173 s_tet[0] = 1.0 / 3.0;
5174 s_tet[1] = 1.0 / 3.0;
5175 s_tet[2] = 1.0 / 3.0;
5176 nod_pt->set_value(0, s_tet[0]);
5177 nod_pt->set_value(1, s_tet[1]);
5178 nod_pt->set_value(2, s_tet[2]);
5179 tet_el_pt->interpolated_x(s_tet, x_tet);
5180 nod_pt->x(0) = x_tet[0];
5181 nod_pt->x(1) = x_tet[1];
5182 nod_pt->x(2) = x_tet[2];
5187 s_tet[0] = 1.0 / 3.0;
5189 s_tet[2] = 1.0 / 3.0;
5190 nod_pt->set_value(0, s_tet[0]);
5191 nod_pt->set_value(1, s_tet[1]);
5192 nod_pt->set_value(2, s_tet[2]);
5193 tet_el_pt->interpolated_x(s_tet, x_tet);
5194 nod_pt->x(0) = x_tet[0];
5195 nod_pt->x(1) = x_tet[1];
5196 nod_pt->x(2) = x_tet[2];
5203 nod_pt->set_value(0, s_tet[0]);
5204 nod_pt->set_value(1, s_tet[1]);
5205 nod_pt->set_value(2, s_tet[2]);
5206 tet_el_pt->interpolated_x(s_tet, x_tet);
5207 nod_pt->x(0) = x_tet[0];
5208 nod_pt->x(1) = x_tet[1];
5209 nod_pt->x(2) = x_tet[2];
5216 FiniteElement* el_pt =
new ELEMENT;
5217 brick_el0_pt = el_pt;
5221 tet_el_pt->node_pt(0), tet_el_pt->node_pt(1), tet_el_pt->node_pt(2));
5224 tet_el_pt->node_pt(0), tet_el_pt->node_pt(2), tet_el_pt->node_pt(3));
5227 tet_el_pt->node_pt(0), tet_el_pt->node_pt(1), tet_el_pt->node_pt(3));
5231 Vector<Vector<unsigned>> tet_edge_node(3);
5232 tet_edge_node[0].resize(2);
5233 tet_edge_node[0][0] = 4;
5234 tet_edge_node[0][1] = 1;
5235 tet_edge_node[1].resize(2);
5236 tet_edge_node[1][0] = 6;
5237 tet_edge_node[1][1] = 3;
5238 tet_edge_node[2].resize(2);
5239 tet_edge_node[2][0] = 5;
5240 tet_edge_node[2][1] = 2;
5243 unsigned central_tet_vertex = 0;
5245 Node* tet_node_pt = 0;
5246 Node* old_node_pt = 0;
5253 tet_node_pt = tet_el_pt->node_pt(central_tet_vertex);
5254 old_node_pt = tet_node_node_pt[tet_node_pt];
5255 if (old_node_pt == 0)
5257 Node* new_node_pt = 0;
5258 if (tet_node_pt->is_on_boundary())
5260 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5264 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5266 tet_node_node_pt[tet_node_pt] = new_node_pt;
5267 Node_pt.push_back(new_node_pt);
5268 Vector<double>
s(3);
5269 Vector<double> s_tet(3);
5270 Vector<double> x_tet(3);
5271 el_pt->local_coordinate_of_node(
j,
s);
5272 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5273 tet_el_pt->interpolated_x(s_tet, x_tet);
5274 new_node_pt->x(0) = x_tet[0];
5275 new_node_pt->x(1) = x_tet[1];
5276 new_node_pt->x(2) = x_tet[2];
5281 el_pt->node_pt(
j) = old_node_pt;
5291 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[0][0]);
5292 old_node_pt = tet_node_node_pt[tet_node_pt];
5293 if (old_node_pt == 0)
5295 Node* new_node_pt = 0;
5296 if (tet_node_pt->is_on_boundary())
5298 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5302 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5304 tet_node_node_pt[tet_node_pt] = new_node_pt;
5305 Node_pt.push_back(new_node_pt);
5306 Vector<double>
s(3);
5307 Vector<double> s_tet(3);
5308 Vector<double> x_tet(3);
5309 el_pt->local_coordinate_of_node(
j,
s);
5310 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5311 tet_el_pt->interpolated_x(s_tet, x_tet);
5312 new_node_pt->x(0) = x_tet[0];
5313 new_node_pt->x(1) = x_tet[1];
5314 new_node_pt->x(2) = x_tet[2];
5319 el_pt->node_pt(
j) = old_node_pt;
5329 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[1][0]);
5330 old_node_pt = tet_node_node_pt[tet_node_pt];
5331 if (old_node_pt == 0)
5333 Node* new_node_pt = 0;
5334 if (tet_node_pt->is_on_boundary())
5336 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5340 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5342 tet_node_node_pt[tet_node_pt] = new_node_pt;
5343 Node_pt.push_back(new_node_pt);
5344 Vector<double>
s(3);
5345 Vector<double> s_tet(3);
5346 Vector<double> x_tet(3);
5347 el_pt->local_coordinate_of_node(
j,
s);
5348 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5349 tet_el_pt->interpolated_x(s_tet, x_tet);
5350 new_node_pt->x(0) = x_tet[0];
5351 new_node_pt->x(1) = x_tet[1];
5352 new_node_pt->x(2) = x_tet[2];
5357 el_pt->node_pt(
j) = old_node_pt;
5367 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[2][0]);
5368 old_node_pt = tet_node_node_pt[tet_node_pt];
5369 if (old_node_pt == 0)
5371 Node* new_node_pt = 0;
5372 if (tet_node_pt->is_on_boundary())
5374 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5378 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5380 tet_node_node_pt[tet_node_pt] = new_node_pt;
5381 Node_pt.push_back(new_node_pt);
5382 Vector<double>
s(3);
5383 Vector<double> s_tet(3);
5384 Vector<double> x_tet(3);
5385 el_pt->local_coordinate_of_node(
j,
s);
5386 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5387 tet_el_pt->interpolated_x(s_tet, x_tet);
5388 new_node_pt->x(0) = x_tet[0];
5389 new_node_pt->x(1) = x_tet[1];
5390 new_node_pt->x(2) = x_tet[2];
5395 el_pt->node_pt(
j) = old_node_pt;
5406 old_node_pt = tet_face_node_pt[
face0];
5407 if (old_node_pt == 0)
5409 Node* new_node_pt = 0;
5410 if (
face0.is_boundary_face())
5412 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5416 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5418 tet_face_node_pt[
face0] = new_node_pt;
5419 Node_pt.push_back(new_node_pt);
5420 Vector<double>
s(3);
5421 Vector<double> s_tet(3);
5422 Vector<double> x_tet(3);
5423 el_pt->local_coordinate_of_node(
j,
s);
5424 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5425 tet_el_pt->interpolated_x(s_tet, x_tet);
5426 new_node_pt->x(0) = x_tet[0];
5427 new_node_pt->x(1) = x_tet[1];
5428 new_node_pt->x(2) = x_tet[2];
5433 el_pt->node_pt(
j) = old_node_pt;
5443 old_node_pt = tet_face_node_pt[
face1];
5444 if (old_node_pt == 0)
5446 Node* new_node_pt = 0;
5447 if (
face1.is_boundary_face())
5449 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5453 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5455 tet_face_node_pt[
face1] = new_node_pt;
5456 Node_pt.push_back(new_node_pt);
5457 Vector<double>
s(3);
5458 Vector<double> s_tet(3);
5459 Vector<double> x_tet(3);
5460 el_pt->local_coordinate_of_node(
j,
s);
5461 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5462 tet_el_pt->interpolated_x(s_tet, x_tet);
5463 new_node_pt->x(0) = x_tet[0];
5464 new_node_pt->x(1) = x_tet[1];
5465 new_node_pt->x(2) = x_tet[2];
5470 el_pt->node_pt(
j) = old_node_pt;
5480 old_node_pt = tet_face_node_pt[
face2];
5481 if (old_node_pt == 0)
5483 Node* new_node_pt = 0;
5484 if (
face2.is_boundary_face())
5486 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5490 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5492 tet_face_node_pt[
face2] = new_node_pt;
5493 Node_pt.push_back(new_node_pt);
5494 Vector<double>
s(3);
5495 Vector<double> s_tet(3);
5496 Vector<double> x_tet(3);
5497 el_pt->local_coordinate_of_node(
j,
s);
5498 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5499 tet_el_pt->interpolated_x(s_tet, x_tet);
5500 new_node_pt->x(0) = x_tet[0];
5501 new_node_pt->x(1) = x_tet[1];
5502 new_node_pt->x(2) = x_tet[2];
5507 el_pt->node_pt(
j) = old_node_pt;
5518 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5519 centroid_node_pt = new_node_pt;
5520 Node_pt.push_back(new_node_pt);
5521 Vector<double>
s(3);
5522 Vector<double> s_tet(3);
5523 Vector<double> x_tet(3);
5524 el_pt->local_coordinate_of_node(
j,
s);
5525 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5526 tet_el_pt->interpolated_x(s_tet, x_tet);
5527 new_node_pt->x(0) = x_tet[0];
5528 new_node_pt->x(1) = x_tet[1];
5529 new_node_pt->x(2) = x_tet[2];
5540 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5541 Node_pt.push_back(new_node_pt);
5542 Vector<double>
s(3);
5543 Vector<double> s_tet(3);
5544 Vector<double> x_tet(3);
5545 el_pt->local_coordinate_of_node(
j,
s);
5546 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5547 tet_el_pt->interpolated_x(s_tet, x_tet);
5548 new_node_pt->x(0) = x_tet[0];
5549 new_node_pt->x(1) = x_tet[1];
5550 new_node_pt->x(2) = x_tet[2];
5559 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
5560 old_node_pt = brick_edge_node_pt[edge];
5561 if (old_node_pt == 0)
5563 Node* new_node_pt = 0;
5564 if (edge.is_boundary_edge())
5566 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5570 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5572 brick_edge_node_pt[edge] = new_node_pt;
5573 Node_pt.push_back(new_node_pt);
5574 Vector<double>
s(3);
5575 Vector<double> s_tet(3);
5576 Vector<double> x_tet(3);
5577 el_pt->local_coordinate_of_node(
j,
s);
5578 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5579 tet_el_pt->interpolated_x(s_tet, x_tet);
5580 new_node_pt->x(0) = x_tet[0];
5581 new_node_pt->x(1) = x_tet[1];
5582 new_node_pt->x(2) = x_tet[2];
5587 el_pt->node_pt(
j) = old_node_pt;
5597 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
5598 old_node_pt = brick_edge_node_pt[edge];
5599 if (old_node_pt == 0)
5601 Node* new_node_pt = 0;
5602 if (edge.is_boundary_edge())
5604 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5608 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5610 brick_edge_node_pt[edge] = new_node_pt;
5611 Node_pt.push_back(new_node_pt);
5612 Vector<double>
s(3);
5613 Vector<double> s_tet(3);
5614 Vector<double> x_tet(3);
5615 el_pt->local_coordinate_of_node(
j,
s);
5616 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5617 tet_el_pt->interpolated_x(s_tet, x_tet);
5618 new_node_pt->x(0) = x_tet[0];
5619 new_node_pt->x(1) = x_tet[1];
5620 new_node_pt->x(2) = x_tet[2];
5625 el_pt->node_pt(
j) = old_node_pt;
5634 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
5635 old_node_pt = brick_edge_node_pt[edge];
5636 if (old_node_pt == 0)
5638 Node* new_node_pt = 0;
5639 if (edge.is_boundary_edge())
5641 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5645 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5647 brick_edge_node_pt[edge] = new_node_pt;
5648 Node_pt.push_back(new_node_pt);
5649 Vector<double>
s(3);
5650 Vector<double> s_tet(3);
5651 Vector<double> x_tet(3);
5652 el_pt->local_coordinate_of_node(
j,
s);
5653 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5654 tet_el_pt->interpolated_x(s_tet, x_tet);
5655 new_node_pt->x(0) = x_tet[0];
5656 new_node_pt->x(1) = x_tet[1];
5657 new_node_pt->x(2) = x_tet[2];
5662 el_pt->node_pt(
j) = old_node_pt;
5671 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
5672 old_node_pt = brick_edge_node_pt[edge];
5673 if (old_node_pt == 0)
5675 Node* new_node_pt = 0;
5676 if (edge.is_boundary_edge())
5678 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5682 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5684 brick_edge_node_pt[edge] = new_node_pt;
5685 Node_pt.push_back(new_node_pt);
5686 Vector<double>
s(3);
5687 Vector<double> s_tet(3);
5688 Vector<double> x_tet(3);
5689 el_pt->local_coordinate_of_node(
j,
s);
5690 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5691 tet_el_pt->interpolated_x(s_tet, x_tet);
5692 new_node_pt->x(0) = x_tet[0];
5693 new_node_pt->x(1) = x_tet[1];
5694 new_node_pt->x(2) = x_tet[2];
5699 el_pt->node_pt(
j) = old_node_pt;
5708 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
5709 old_node_pt = brick_edge_node_pt[edge];
5710 if (old_node_pt == 0)
5712 Node* new_node_pt = 0;
5713 if (edge.is_boundary_edge())
5715 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5719 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5721 brick_edge_node_pt[edge] = new_node_pt;
5722 Node_pt.push_back(new_node_pt);
5723 Vector<double>
s(3);
5724 Vector<double> s_tet(3);
5725 Vector<double> x_tet(3);
5726 el_pt->local_coordinate_of_node(
j,
s);
5727 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5728 tet_el_pt->interpolated_x(s_tet, x_tet);
5729 new_node_pt->x(0) = x_tet[0];
5730 new_node_pt->x(1) = x_tet[1];
5731 new_node_pt->x(2) = x_tet[2];
5736 el_pt->node_pt(
j) = old_node_pt;
5746 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
5747 old_node_pt = brick_edge_node_pt[edge];
5748 if (old_node_pt == 0)
5750 Node* new_node_pt = 0;
5751 if (edge.is_boundary_edge())
5753 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5757 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5759 brick_edge_node_pt[edge] = new_node_pt;
5760 Node_pt.push_back(new_node_pt);
5761 Vector<double>
s(3);
5762 Vector<double> s_tet(3);
5763 Vector<double> x_tet(3);
5764 el_pt->local_coordinate_of_node(
j,
s);
5765 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5766 tet_el_pt->interpolated_x(s_tet, x_tet);
5767 new_node_pt->x(0) = x_tet[0];
5768 new_node_pt->x(1) = x_tet[1];
5769 new_node_pt->x(2) = x_tet[2];
5774 el_pt->node_pt(
j) = old_node_pt;
5783 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
5784 old_node_pt = brick_edge_node_pt[edge];
5785 if (old_node_pt == 0)
5787 Node* new_node_pt = 0;
5788 if (edge.is_boundary_edge())
5790 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5794 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5796 brick_edge_node_pt[edge] = new_node_pt;
5797 Node_pt.push_back(new_node_pt);
5798 Vector<double>
s(3);
5799 Vector<double> s_tet(3);
5800 Vector<double> x_tet(3);
5801 el_pt->local_coordinate_of_node(
j,
s);
5802 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5803 tet_el_pt->interpolated_x(s_tet, x_tet);
5804 new_node_pt->x(0) = x_tet[0];
5805 new_node_pt->x(1) = x_tet[1];
5806 new_node_pt->x(2) = x_tet[2];
5811 el_pt->node_pt(
j) = old_node_pt;
5821 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
5822 old_node_pt = brick_edge_node_pt[edge];
5823 if (old_node_pt == 0)
5825 Node* new_node_pt = 0;
5826 if (edge.is_boundary_edge())
5828 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5832 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5834 brick_edge_node_pt[edge] = new_node_pt;
5835 Node_pt.push_back(new_node_pt);
5836 Vector<double>
s(3);
5837 Vector<double> s_tet(3);
5838 Vector<double> x_tet(3);
5839 el_pt->local_coordinate_of_node(
j,
s);
5840 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5841 tet_el_pt->interpolated_x(s_tet, x_tet);
5842 new_node_pt->x(0) = x_tet[0];
5843 new_node_pt->x(1) = x_tet[1];
5844 new_node_pt->x(2) = x_tet[2];
5849 el_pt->node_pt(
j) = old_node_pt;
5858 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
5859 old_node_pt = brick_edge_node_pt[edge];
5860 if (old_node_pt == 0)
5862 Node* new_node_pt = 0;
5863 if (edge.is_boundary_edge())
5865 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5869 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5871 brick_edge_node_pt[edge] = new_node_pt;
5872 Node_pt.push_back(new_node_pt);
5873 Vector<double>
s(3);
5874 Vector<double> s_tet(3);
5875 Vector<double> x_tet(3);
5876 el_pt->local_coordinate_of_node(
j,
s);
5877 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5878 tet_el_pt->interpolated_x(s_tet, x_tet);
5879 new_node_pt->x(0) = x_tet[0];
5880 new_node_pt->x(1) = x_tet[1];
5881 new_node_pt->x(2) = x_tet[2];
5886 el_pt->node_pt(
j) = old_node_pt;
5896 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
5897 old_node_pt = brick_edge_node_pt[edge];
5898 if (old_node_pt == 0)
5900 Node* new_node_pt = 0;
5901 if (edge.is_boundary_edge())
5903 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5907 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5909 brick_edge_node_pt[edge] = new_node_pt;
5910 Node_pt.push_back(new_node_pt);
5911 Vector<double>
s(3);
5912 Vector<double> s_tet(3);
5913 Vector<double> x_tet(3);
5914 el_pt->local_coordinate_of_node(
j,
s);
5915 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5916 tet_el_pt->interpolated_x(s_tet, x_tet);
5917 new_node_pt->x(0) = x_tet[0];
5918 new_node_pt->x(1) = x_tet[1];
5919 new_node_pt->x(2) = x_tet[2];
5924 el_pt->node_pt(
j) = old_node_pt;
5934 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
5935 old_node_pt = brick_edge_node_pt[edge];
5936 if (old_node_pt == 0)
5938 Node* new_node_pt = 0;
5939 if (edge.is_boundary_edge())
5941 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5945 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5947 brick_edge_node_pt[edge] = new_node_pt;
5948 Node_pt.push_back(new_node_pt);
5949 Vector<double>
s(3);
5950 Vector<double> s_tet(3);
5951 Vector<double> x_tet(3);
5952 el_pt->local_coordinate_of_node(
j,
s);
5953 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5954 tet_el_pt->interpolated_x(s_tet, x_tet);
5955 new_node_pt->x(0) = x_tet[0];
5956 new_node_pt->x(1) = x_tet[1];
5957 new_node_pt->x(2) = x_tet[2];
5962 el_pt->node_pt(
j) = old_node_pt;
5972 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
5973 old_node_pt = brick_edge_node_pt[edge];
5974 if (old_node_pt == 0)
5976 Node* new_node_pt = 0;
5977 if (edge.is_boundary_edge())
5979 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
5983 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
5985 brick_edge_node_pt[edge] = new_node_pt;
5986 Node_pt.push_back(new_node_pt);
5987 Vector<double>
s(3);
5988 Vector<double> s_tet(3);
5989 Vector<double> x_tet(3);
5990 el_pt->local_coordinate_of_node(
j,
s);
5991 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
5992 tet_el_pt->interpolated_x(s_tet, x_tet);
5993 new_node_pt->x(0) = x_tet[0];
5994 new_node_pt->x(1) = x_tet[1];
5995 new_node_pt->x(2) = x_tet[2];
6000 el_pt->node_pt(
j) = old_node_pt;
6011 TFace face(tet_el_pt->node_pt(central_tet_vertex),
6012 tet_el_pt->node_pt(tet_edge_node[0][0]),
6013 tet_el_pt->node_pt(tet_edge_node[2][0]));
6015 old_node_pt = tet_face_node_pt[face];
6016 if (old_node_pt == 0)
6018 Node* new_node_pt = 0;
6019 if (face.is_boundary_face())
6021 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6025 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6027 tet_face_node_pt[face] = new_node_pt;
6028 Node_pt.push_back(new_node_pt);
6029 Vector<double>
s(3);
6030 Vector<double> s_tet(3);
6031 Vector<double> x_tet(3);
6032 el_pt->local_coordinate_of_node(
j,
s);
6033 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
6034 tet_el_pt->interpolated_x(s_tet, x_tet);
6035 new_node_pt->x(0) = x_tet[0];
6036 new_node_pt->x(1) = x_tet[1];
6037 new_node_pt->x(2) = x_tet[2];
6042 el_pt->node_pt(
j) = old_node_pt;
6053 TFace face(tet_el_pt->node_pt(central_tet_vertex),
6054 tet_el_pt->node_pt(tet_edge_node[1][0]),
6055 tet_el_pt->node_pt(tet_edge_node[2][0]));
6057 old_node_pt = tet_face_node_pt[face];
6058 if (old_node_pt == 0)
6060 Node* new_node_pt = 0;
6061 if (face.is_boundary_face())
6063 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6067 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6069 tet_face_node_pt[face] = new_node_pt;
6070 Node_pt.push_back(new_node_pt);
6071 Vector<double>
s(3);
6072 Vector<double> s_tet(3);
6073 Vector<double> x_tet(3);
6074 el_pt->local_coordinate_of_node(
j,
s);
6075 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
6076 tet_el_pt->interpolated_x(s_tet, x_tet);
6077 new_node_pt->x(0) = x_tet[0];
6078 new_node_pt->x(1) = x_tet[1];
6079 new_node_pt->x(2) = x_tet[2];
6084 el_pt->node_pt(
j) = old_node_pt;
6095 TFace face(tet_el_pt->node_pt(central_tet_vertex),
6096 tet_el_pt->node_pt(tet_edge_node[0][0]),
6097 tet_el_pt->node_pt(tet_edge_node[1][0]));
6099 old_node_pt = tet_face_node_pt[face];
6100 if (old_node_pt == 0)
6102 Node* new_node_pt = 0;
6103 if (face.is_boundary_face())
6105 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6109 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6111 tet_face_node_pt[face] = new_node_pt;
6112 Node_pt.push_back(new_node_pt);
6113 Vector<double>
s(3);
6114 Vector<double> s_tet(3);
6115 Vector<double> x_tet(3);
6116 el_pt->local_coordinate_of_node(
j,
s);
6117 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
6118 tet_el_pt->interpolated_x(s_tet, x_tet);
6119 new_node_pt->x(0) = x_tet[0];
6120 new_node_pt->x(1) = x_tet[1];
6121 new_node_pt->x(2) = x_tet[2];
6126 el_pt->node_pt(
j) = old_node_pt;
6134 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6135 Node_pt.push_back(new_node_pt);
6136 Vector<double>
s(3);
6137 Vector<double> s_tet(3);
6138 Vector<double> x_tet(3);
6139 el_pt->local_coordinate_of_node(
j,
s);
6140 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
6141 top_mid_face_node0_pt = new_node_pt;
6142 tet_el_pt->interpolated_x(s_tet, x_tet);
6143 new_node_pt->x(0) = x_tet[0];
6144 new_node_pt->x(1) = x_tet[1];
6145 new_node_pt->x(2) = x_tet[2];
6152 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6153 Node_pt.push_back(new_node_pt);
6154 Vector<double>
s(3);
6155 Vector<double> s_tet(3);
6156 Vector<double> x_tet(3);
6157 el_pt->local_coordinate_of_node(
j,
s);
6158 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
6159 right_mid_face_node0_pt = new_node_pt;
6160 tet_el_pt->interpolated_x(s_tet, x_tet);
6161 new_node_pt->x(0) = x_tet[0];
6162 new_node_pt->x(1) = x_tet[1];
6163 new_node_pt->x(2) = x_tet[2];
6170 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6171 Node_pt.push_back(new_node_pt);
6172 Vector<double>
s(3);
6173 Vector<double> s_tet(3);
6174 Vector<double> x_tet(3);
6175 el_pt->local_coordinate_of_node(
j,
s);
6176 dummy_q_el_pt[0]->interpolated_s_tet(
s, s_tet);
6177 back_mid_face_node0_pt = new_node_pt;
6178 tet_el_pt->interpolated_x(s_tet, x_tet);
6179 new_node_pt->x(0) = x_tet[0];
6180 new_node_pt->x(1) = x_tet[1];
6181 new_node_pt->x(2) = x_tet[2];
6190 for (
unsigned j = 0;
j < 8;
j++)
6192 Node* nod_pt = dummy_q_el_pt[1]->node_pt(
j);
6193 Vector<double> s_tet(3);
6194 Vector<double> x_tet(3);
6198 tet_el_pt->local_coordinate_of_node(1, s_tet);
6199 nod_pt->set_value(0, s_tet[0]);
6200 nod_pt->set_value(1, s_tet[1]);
6201 nod_pt->set_value(2, s_tet[2]);
6202 tet_el_pt->interpolated_x(s_tet, x_tet);
6203 nod_pt->x(0) = x_tet[0];
6204 nod_pt->x(1) = x_tet[1];
6205 nod_pt->x(2) = x_tet[2];
6208 tet_el_pt->local_coordinate_of_node(9, s_tet);
6209 nod_pt->set_value(0, s_tet[0]);
6210 nod_pt->set_value(1, s_tet[1]);
6211 nod_pt->set_value(2, s_tet[2]);
6212 tet_el_pt->interpolated_x(s_tet, x_tet);
6213 nod_pt->x(0) = x_tet[0];
6214 nod_pt->x(1) = x_tet[1];
6215 nod_pt->x(2) = x_tet[2];
6218 tet_el_pt->local_coordinate_of_node(4, s_tet);
6219 nod_pt->set_value(0, s_tet[0]);
6220 nod_pt->set_value(1, s_tet[1]);
6221 nod_pt->set_value(2, s_tet[2]);
6222 tet_el_pt->interpolated_x(s_tet, x_tet);
6223 nod_pt->x(0) = x_tet[0];
6224 nod_pt->x(1) = x_tet[1];
6225 nod_pt->x(2) = x_tet[2];
6230 s_tet[0] = 1.0 / 3.0;
6231 s_tet[1] = 1.0 / 3.0;
6233 nod_pt->set_value(0, s_tet[0]);
6234 nod_pt->set_value(1, s_tet[1]);
6235 nod_pt->set_value(2, s_tet[2]);
6236 tet_el_pt->interpolated_x(s_tet, x_tet);
6237 nod_pt->x(0) = x_tet[0];
6238 nod_pt->x(1) = x_tet[1];
6239 nod_pt->x(2) = x_tet[2];
6242 tet_el_pt->local_coordinate_of_node(7, s_tet);
6243 nod_pt->set_value(0, s_tet[0]);
6244 nod_pt->set_value(1, s_tet[1]);
6245 nod_pt->set_value(2, s_tet[2]);
6246 tet_el_pt->interpolated_x(s_tet, x_tet);
6247 nod_pt->x(0) = x_tet[0];
6248 nod_pt->x(1) = x_tet[1];
6249 nod_pt->x(2) = x_tet[2];
6255 s_tet[1] = 1.0 / 3.0;
6256 s_tet[2] = 1.0 / 3.0;
6257 nod_pt->set_value(0, s_tet[0]);
6258 nod_pt->set_value(1, s_tet[1]);
6259 nod_pt->set_value(2, s_tet[2]);
6260 tet_el_pt->interpolated_x(s_tet, x_tet);
6261 nod_pt->x(0) = x_tet[0];
6262 nod_pt->x(1) = x_tet[1];
6263 nod_pt->x(2) = x_tet[2];
6268 s_tet[0] = 1.0 / 3.0;
6269 s_tet[1] = 1.0 / 3.0;
6270 s_tet[2] = 1.0 / 3.0;
6271 nod_pt->set_value(0, s_tet[0]);
6272 nod_pt->set_value(1, s_tet[1]);
6273 nod_pt->set_value(2, s_tet[2]);
6274 tet_el_pt->interpolated_x(s_tet, x_tet);
6275 nod_pt->x(0) = x_tet[0];
6276 nod_pt->x(1) = x_tet[1];
6277 nod_pt->x(2) = x_tet[2];
6284 nod_pt->set_value(0, s_tet[0]);
6285 nod_pt->set_value(1, s_tet[1]);
6286 nod_pt->set_value(2, s_tet[2]);
6287 tet_el_pt->interpolated_x(s_tet, x_tet);
6288 nod_pt->x(0) = x_tet[0];
6289 nod_pt->x(1) = x_tet[1];
6290 nod_pt->x(2) = x_tet[2];
6297 FiniteElement* el_pt =
new ELEMENT;
6298 brick_el1_pt = el_pt;
6302 tet_el_pt->node_pt(1), tet_el_pt->node_pt(3), tet_el_pt->node_pt(2));
6305 tet_el_pt->node_pt(1), tet_el_pt->node_pt(0), tet_el_pt->node_pt(2));
6308 tet_el_pt->node_pt(1), tet_el_pt->node_pt(0), tet_el_pt->node_pt(3));
6311 Vector<Vector<unsigned>> tet_edge_node(3);
6312 tet_edge_node[0].resize(2);
6313 tet_edge_node[0][0] = 9;
6314 tet_edge_node[0][1] = 3;
6315 tet_edge_node[1].resize(2);
6316 tet_edge_node[1][0] = 4;
6317 tet_edge_node[1][1] = 0;
6318 tet_edge_node[2].resize(2);
6319 tet_edge_node[2][0] = 7;
6320 tet_edge_node[2][1] = 2;
6323 unsigned central_tet_vertex = 1;
6325 Node* tet_node_pt = 0;
6326 Node* old_node_pt = 0;
6333 tet_node_pt = tet_el_pt->node_pt(central_tet_vertex);
6334 old_node_pt = tet_node_node_pt[tet_node_pt];
6335 if (old_node_pt == 0)
6337 Node* new_node_pt = 0;
6338 if (tet_node_pt->is_on_boundary())
6340 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6344 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6346 tet_node_node_pt[tet_node_pt] = new_node_pt;
6347 Node_pt.push_back(new_node_pt);
6348 Vector<double>
s(3);
6349 Vector<double> s_tet(3);
6350 Vector<double> x_tet(3);
6351 el_pt->local_coordinate_of_node(
j,
s);
6352 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6353 tet_el_pt->interpolated_x(s_tet, x_tet);
6354 new_node_pt->x(0) = x_tet[0];
6355 new_node_pt->x(1) = x_tet[1];
6356 new_node_pt->x(2) = x_tet[2];
6361 el_pt->node_pt(
j) = old_node_pt;
6371 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[0][0]);
6372 old_node_pt = tet_node_node_pt[tet_node_pt];
6373 if (old_node_pt == 0)
6375 Node* new_node_pt = 0;
6376 if (tet_node_pt->is_on_boundary())
6378 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6382 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6384 tet_node_node_pt[tet_node_pt] = new_node_pt;
6385 Node_pt.push_back(new_node_pt);
6386 Vector<double>
s(3);
6387 Vector<double> s_tet(3);
6388 Vector<double> x_tet(3);
6389 el_pt->local_coordinate_of_node(
j,
s);
6390 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6391 tet_el_pt->interpolated_x(s_tet, x_tet);
6392 new_node_pt->x(0) = x_tet[0];
6393 new_node_pt->x(1) = x_tet[1];
6394 new_node_pt->x(2) = x_tet[2];
6399 el_pt->node_pt(
j) = old_node_pt;
6409 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[1][0]);
6410 old_node_pt = tet_node_node_pt[tet_node_pt];
6411 if (old_node_pt == 0)
6413 Node* new_node_pt = 0;
6414 if (tet_node_pt->is_on_boundary())
6416 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6420 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6422 tet_node_node_pt[tet_node_pt] = new_node_pt;
6423 Node_pt.push_back(new_node_pt);
6424 Vector<double>
s(3);
6425 Vector<double> s_tet(3);
6426 Vector<double> x_tet(3);
6427 el_pt->local_coordinate_of_node(
j,
s);
6428 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6429 tet_el_pt->interpolated_x(s_tet, x_tet);
6430 new_node_pt->x(0) = x_tet[0];
6431 new_node_pt->x(1) = x_tet[1];
6432 new_node_pt->x(2) = x_tet[2];
6437 el_pt->node_pt(
j) = old_node_pt;
6447 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[2][0]);
6448 old_node_pt = tet_node_node_pt[tet_node_pt];
6449 if (old_node_pt == 0)
6451 Node* new_node_pt = 0;
6452 if (tet_node_pt->is_on_boundary())
6454 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6458 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6460 tet_node_node_pt[tet_node_pt] = new_node_pt;
6461 Node_pt.push_back(new_node_pt);
6462 Vector<double>
s(3);
6463 Vector<double> s_tet(3);
6464 Vector<double> x_tet(3);
6465 el_pt->local_coordinate_of_node(
j,
s);
6466 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6467 tet_el_pt->interpolated_x(s_tet, x_tet);
6468 new_node_pt->x(0) = x_tet[0];
6469 new_node_pt->x(1) = x_tet[1];
6470 new_node_pt->x(2) = x_tet[2];
6475 el_pt->node_pt(
j) = old_node_pt;
6485 old_node_pt = tet_face_node_pt[
face0];
6486 if (old_node_pt == 0)
6488 Node* new_node_pt = 0;
6489 if (
face0.is_boundary_face())
6491 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6495 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6497 tet_face_node_pt[
face0] = new_node_pt;
6498 Node_pt.push_back(new_node_pt);
6499 Vector<double>
s(3);
6500 Vector<double> s_tet(3);
6501 Vector<double> x_tet(3);
6502 el_pt->local_coordinate_of_node(
j,
s);
6503 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6504 tet_el_pt->interpolated_x(s_tet, x_tet);
6505 new_node_pt->x(0) = x_tet[0];
6506 new_node_pt->x(1) = x_tet[1];
6507 new_node_pt->x(2) = x_tet[2];
6512 el_pt->node_pt(
j) = old_node_pt;
6521 old_node_pt = tet_face_node_pt[
face1];
6522 if (old_node_pt == 0)
6524 Node* new_node_pt = 0;
6525 if (
face1.is_boundary_face())
6527 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6531 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6533 tet_face_node_pt[
face1] = new_node_pt;
6534 Node_pt.push_back(new_node_pt);
6535 Vector<double>
s(3);
6536 Vector<double> s_tet(3);
6537 Vector<double> x_tet(3);
6538 el_pt->local_coordinate_of_node(
j,
s);
6539 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6540 tet_el_pt->interpolated_x(s_tet, x_tet);
6541 new_node_pt->x(0) = x_tet[0];
6542 new_node_pt->x(1) = x_tet[1];
6543 new_node_pt->x(2) = x_tet[2];
6548 el_pt->node_pt(
j) = old_node_pt;
6557 old_node_pt = tet_face_node_pt[
face2];
6558 if (old_node_pt == 0)
6560 Node* new_node_pt = 0;
6561 if (
face2.is_boundary_face())
6563 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6567 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6569 tet_face_node_pt[
face2] = new_node_pt;
6570 Node_pt.push_back(new_node_pt);
6571 Vector<double>
s(3);
6572 Vector<double> s_tet(3);
6573 Vector<double> x_tet(3);
6574 el_pt->local_coordinate_of_node(
j,
s);
6575 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6576 tet_el_pt->interpolated_x(s_tet, x_tet);
6577 new_node_pt->x(0) = x_tet[0];
6578 new_node_pt->x(1) = x_tet[1];
6579 new_node_pt->x(2) = x_tet[2];
6584 el_pt->node_pt(
j) = old_node_pt;
6595 el_pt->node_pt(
j) = centroid_node_pt;
6605 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6606 Node_pt.push_back(new_node_pt);
6607 Vector<double>
s(3);
6608 Vector<double> s_tet(3);
6609 Vector<double> x_tet(3);
6610 el_pt->local_coordinate_of_node(
j,
s);
6611 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6612 tet_el_pt->interpolated_x(s_tet, x_tet);
6613 new_node_pt->x(0) = x_tet[0];
6614 new_node_pt->x(1) = x_tet[1];
6615 new_node_pt->x(2) = x_tet[2];
6625 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
6626 old_node_pt = brick_edge_node_pt[edge];
6627 if (old_node_pt == 0)
6629 Node* new_node_pt = 0;
6630 if (edge.is_boundary_edge())
6632 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6636 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6638 brick_edge_node_pt[edge] = new_node_pt;
6639 Node_pt.push_back(new_node_pt);
6640 Vector<double>
s(3);
6641 Vector<double> s_tet(3);
6642 Vector<double> x_tet(3);
6643 el_pt->local_coordinate_of_node(
j,
s);
6644 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6645 tet_el_pt->interpolated_x(s_tet, x_tet);
6646 new_node_pt->x(0) = x_tet[0];
6647 new_node_pt->x(1) = x_tet[1];
6648 new_node_pt->x(2) = x_tet[2];
6653 el_pt->node_pt(
j) = old_node_pt;
6663 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
6664 old_node_pt = brick_edge_node_pt[edge];
6665 if (old_node_pt == 0)
6667 Node* new_node_pt = 0;
6668 if (edge.is_boundary_edge())
6670 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6674 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6676 brick_edge_node_pt[edge] = new_node_pt;
6677 Node_pt.push_back(new_node_pt);
6678 Vector<double>
s(3);
6679 Vector<double> s_tet(3);
6680 Vector<double> x_tet(3);
6681 el_pt->local_coordinate_of_node(
j,
s);
6682 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6683 tet_el_pt->interpolated_x(s_tet, x_tet);
6684 new_node_pt->x(0) = x_tet[0];
6685 new_node_pt->x(1) = x_tet[1];
6686 new_node_pt->x(2) = x_tet[2];
6691 el_pt->node_pt(
j) = old_node_pt;
6701 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
6702 old_node_pt = brick_edge_node_pt[edge];
6703 if (old_node_pt == 0)
6705 Node* new_node_pt = 0;
6706 if (edge.is_boundary_edge())
6708 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6712 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6714 brick_edge_node_pt[edge] = new_node_pt;
6715 Node_pt.push_back(new_node_pt);
6716 Vector<double>
s(3);
6717 Vector<double> s_tet(3);
6718 Vector<double> x_tet(3);
6719 el_pt->local_coordinate_of_node(
j,
s);
6720 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6721 tet_el_pt->interpolated_x(s_tet, x_tet);
6722 new_node_pt->x(0) = x_tet[0];
6723 new_node_pt->x(1) = x_tet[1];
6724 new_node_pt->x(2) = x_tet[2];
6729 el_pt->node_pt(
j) = old_node_pt;
6738 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
6739 old_node_pt = brick_edge_node_pt[edge];
6740 if (old_node_pt == 0)
6742 Node* new_node_pt = 0;
6743 if (edge.is_boundary_edge())
6745 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6749 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6751 brick_edge_node_pt[edge] = new_node_pt;
6752 Node_pt.push_back(new_node_pt);
6753 Vector<double>
s(3);
6754 Vector<double> s_tet(3);
6755 Vector<double> x_tet(3);
6756 el_pt->local_coordinate_of_node(
j,
s);
6757 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6758 tet_el_pt->interpolated_x(s_tet, x_tet);
6759 new_node_pt->x(0) = x_tet[0];
6760 new_node_pt->x(1) = x_tet[1];
6761 new_node_pt->x(2) = x_tet[2];
6766 el_pt->node_pt(
j) = old_node_pt;
6775 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
6776 old_node_pt = brick_edge_node_pt[edge];
6777 if (old_node_pt == 0)
6779 Node* new_node_pt = 0;
6780 if (edge.is_boundary_edge())
6782 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6786 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6788 brick_edge_node_pt[edge] = new_node_pt;
6789 Node_pt.push_back(new_node_pt);
6790 Vector<double>
s(3);
6791 Vector<double> s_tet(3);
6792 Vector<double> x_tet(3);
6793 el_pt->local_coordinate_of_node(
j,
s);
6794 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6795 tet_el_pt->interpolated_x(s_tet, x_tet);
6796 new_node_pt->x(0) = x_tet[0];
6797 new_node_pt->x(1) = x_tet[1];
6798 new_node_pt->x(2) = x_tet[2];
6803 el_pt->node_pt(
j) = old_node_pt;
6813 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
6814 old_node_pt = brick_edge_node_pt[edge];
6815 if (old_node_pt == 0)
6817 Node* new_node_pt = 0;
6818 if (edge.is_boundary_edge())
6820 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6824 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6826 brick_edge_node_pt[edge] = new_node_pt;
6827 Node_pt.push_back(new_node_pt);
6828 Vector<double>
s(3);
6829 Vector<double> s_tet(3);
6830 Vector<double> x_tet(3);
6831 el_pt->local_coordinate_of_node(
j,
s);
6832 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6833 tet_el_pt->interpolated_x(s_tet, x_tet);
6834 new_node_pt->x(0) = x_tet[0];
6835 new_node_pt->x(1) = x_tet[1];
6836 new_node_pt->x(2) = x_tet[2];
6841 el_pt->node_pt(
j) = old_node_pt;
6850 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
6851 old_node_pt = brick_edge_node_pt[edge];
6852 if (old_node_pt == 0)
6854 Node* new_node_pt = 0;
6855 if (edge.is_boundary_edge())
6857 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6861 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6863 brick_edge_node_pt[edge] = new_node_pt;
6864 Node_pt.push_back(new_node_pt);
6865 Vector<double>
s(3);
6866 Vector<double> s_tet(3);
6867 Vector<double> x_tet(3);
6868 el_pt->local_coordinate_of_node(
j,
s);
6869 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6870 tet_el_pt->interpolated_x(s_tet, x_tet);
6871 new_node_pt->x(0) = x_tet[0];
6872 new_node_pt->x(1) = x_tet[1];
6873 new_node_pt->x(2) = x_tet[2];
6878 el_pt->node_pt(
j) = old_node_pt;
6888 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
6889 old_node_pt = brick_edge_node_pt[edge];
6890 if (old_node_pt == 0)
6892 Node* new_node_pt = 0;
6893 if (edge.is_boundary_edge())
6895 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6899 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6901 brick_edge_node_pt[edge] = new_node_pt;
6902 Node_pt.push_back(new_node_pt);
6903 Vector<double>
s(3);
6904 Vector<double> s_tet(3);
6905 Vector<double> x_tet(3);
6906 el_pt->local_coordinate_of_node(
j,
s);
6907 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6908 tet_el_pt->interpolated_x(s_tet, x_tet);
6909 new_node_pt->x(0) = x_tet[0];
6910 new_node_pt->x(1) = x_tet[1];
6911 new_node_pt->x(2) = x_tet[2];
6916 el_pt->node_pt(
j) = old_node_pt;
6925 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
6926 old_node_pt = brick_edge_node_pt[edge];
6927 if (old_node_pt == 0)
6929 Node* new_node_pt = 0;
6930 if (edge.is_boundary_edge())
6932 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6936 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6938 brick_edge_node_pt[edge] = new_node_pt;
6939 Node_pt.push_back(new_node_pt);
6940 Vector<double>
s(3);
6941 Vector<double> s_tet(3);
6942 Vector<double> x_tet(3);
6943 el_pt->local_coordinate_of_node(
j,
s);
6944 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6945 tet_el_pt->interpolated_x(s_tet, x_tet);
6946 new_node_pt->x(0) = x_tet[0];
6947 new_node_pt->x(1) = x_tet[1];
6948 new_node_pt->x(2) = x_tet[2];
6953 el_pt->node_pt(
j) = old_node_pt;
6963 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
6964 old_node_pt = brick_edge_node_pt[edge];
6965 if (old_node_pt == 0)
6967 Node* new_node_pt = 0;
6968 if (edge.is_boundary_edge())
6970 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
6974 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
6976 brick_edge_node_pt[edge] = new_node_pt;
6977 Node_pt.push_back(new_node_pt);
6978 Vector<double>
s(3);
6979 Vector<double> s_tet(3);
6980 Vector<double> x_tet(3);
6981 el_pt->local_coordinate_of_node(
j,
s);
6982 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
6983 tet_el_pt->interpolated_x(s_tet, x_tet);
6984 new_node_pt->x(0) = x_tet[0];
6985 new_node_pt->x(1) = x_tet[1];
6986 new_node_pt->x(2) = x_tet[2];
6991 el_pt->node_pt(
j) = old_node_pt;
7001 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
7002 old_node_pt = brick_edge_node_pt[edge];
7003 if (old_node_pt == 0)
7005 Node* new_node_pt = 0;
7006 if (edge.is_boundary_edge())
7008 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7012 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7014 brick_edge_node_pt[edge] = new_node_pt;
7015 Node_pt.push_back(new_node_pt);
7016 Vector<double>
s(3);
7017 Vector<double> s_tet(3);
7018 Vector<double> x_tet(3);
7019 el_pt->local_coordinate_of_node(
j,
s);
7020 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
7021 tet_el_pt->interpolated_x(s_tet, x_tet);
7022 new_node_pt->x(0) = x_tet[0];
7023 new_node_pt->x(1) = x_tet[1];
7024 new_node_pt->x(2) = x_tet[2];
7029 el_pt->node_pt(
j) = old_node_pt;
7039 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
7040 old_node_pt = brick_edge_node_pt[edge];
7041 if (old_node_pt == 0)
7043 Node* new_node_pt = 0;
7044 if (edge.is_boundary_edge())
7046 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7050 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7052 brick_edge_node_pt[edge] = new_node_pt;
7053 Node_pt.push_back(new_node_pt);
7054 Vector<double>
s(3);
7055 Vector<double> s_tet(3);
7056 Vector<double> x_tet(3);
7057 el_pt->local_coordinate_of_node(
j,
s);
7058 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
7059 tet_el_pt->interpolated_x(s_tet, x_tet);
7060 new_node_pt->x(0) = x_tet[0];
7061 new_node_pt->x(1) = x_tet[1];
7062 new_node_pt->x(2) = x_tet[2];
7067 el_pt->node_pt(
j) = old_node_pt;
7078 TFace face(tet_el_pt->node_pt(central_tet_vertex),
7079 tet_el_pt->node_pt(tet_edge_node[0][0]),
7080 tet_el_pt->node_pt(tet_edge_node[2][0]));
7082 old_node_pt = tet_face_node_pt[face];
7083 if (old_node_pt == 0)
7085 Node* new_node_pt = 0;
7086 if (face.is_boundary_face())
7088 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7092 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7094 tet_face_node_pt[face] = new_node_pt;
7095 Node_pt.push_back(new_node_pt);
7096 Vector<double>
s(3);
7097 Vector<double> s_tet(3);
7098 Vector<double> x_tet(3);
7099 el_pt->local_coordinate_of_node(
j,
s);
7100 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
7101 tet_el_pt->interpolated_x(s_tet, x_tet);
7102 new_node_pt->x(0) = x_tet[0];
7103 new_node_pt->x(1) = x_tet[1];
7104 new_node_pt->x(2) = x_tet[2];
7109 el_pt->node_pt(
j) = old_node_pt;
7120 TFace face(tet_el_pt->node_pt(central_tet_vertex),
7121 tet_el_pt->node_pt(tet_edge_node[1][0]),
7122 tet_el_pt->node_pt(tet_edge_node[2][0]));
7124 old_node_pt = tet_face_node_pt[face];
7125 if (old_node_pt == 0)
7127 Node* new_node_pt = 0;
7128 if (face.is_boundary_face())
7130 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7134 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7136 tet_face_node_pt[face] = new_node_pt;
7137 Node_pt.push_back(new_node_pt);
7138 Vector<double>
s(3);
7139 Vector<double> s_tet(3);
7140 Vector<double> x_tet(3);
7141 el_pt->local_coordinate_of_node(
j,
s);
7142 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
7143 tet_el_pt->interpolated_x(s_tet, x_tet);
7144 new_node_pt->x(0) = x_tet[0];
7145 new_node_pt->x(1) = x_tet[1];
7146 new_node_pt->x(2) = x_tet[2];
7151 el_pt->node_pt(
j) = old_node_pt;
7162 TFace face(tet_el_pt->node_pt(central_tet_vertex),
7163 tet_el_pt->node_pt(tet_edge_node[0][0]),
7164 tet_el_pt->node_pt(tet_edge_node[1][0]));
7166 old_node_pt = tet_face_node_pt[face];
7167 if (old_node_pt == 0)
7169 Node* new_node_pt = 0;
7170 if (face.is_boundary_face())
7172 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7176 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7178 tet_face_node_pt[face] = new_node_pt;
7179 Node_pt.push_back(new_node_pt);
7180 Vector<double>
s(3);
7181 Vector<double> s_tet(3);
7182 Vector<double> x_tet(3);
7183 el_pt->local_coordinate_of_node(
j,
s);
7184 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
7185 tet_el_pt->interpolated_x(s_tet, x_tet);
7186 new_node_pt->x(0) = x_tet[0];
7187 new_node_pt->x(1) = x_tet[1];
7188 new_node_pt->x(2) = x_tet[2];
7193 el_pt->node_pt(
j) = old_node_pt;
7201 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7202 Node_pt.push_back(new_node_pt);
7203 Vector<double>
s(3);
7204 Vector<double> s_tet(3);
7205 Vector<double> x_tet(3);
7206 el_pt->local_coordinate_of_node(
j,
s);
7207 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
7208 top_mid_face_node1_pt = new_node_pt;
7209 tet_el_pt->interpolated_x(s_tet, x_tet);
7210 new_node_pt->x(0) = x_tet[0];
7211 new_node_pt->x(1) = x_tet[1];
7212 new_node_pt->x(2) = x_tet[2];
7219 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7220 Node_pt.push_back(new_node_pt);
7221 Vector<double>
s(3);
7222 Vector<double> s_tet(3);
7223 Vector<double> x_tet(3);
7224 el_pt->local_coordinate_of_node(
j,
s);
7225 dummy_q_el_pt[1]->interpolated_s_tet(
s, s_tet);
7226 right_mid_face_node1_pt = new_node_pt;
7227 tet_el_pt->interpolated_x(s_tet, x_tet);
7228 new_node_pt->x(0) = x_tet[0];
7229 new_node_pt->x(1) = x_tet[1];
7230 new_node_pt->x(2) = x_tet[2];
7239 el_pt->node_pt(
j) = right_mid_face_node0_pt;
7248 for (
unsigned j = 0;
j < 8;
j++)
7250 Node* nod_pt = dummy_q_el_pt[2]->node_pt(
j);
7251 Vector<double> s_tet(3);
7252 Vector<double> x_tet(3);
7256 tet_el_pt->local_coordinate_of_node(3, s_tet);
7257 nod_pt->set_value(0, s_tet[0]);
7258 nod_pt->set_value(1, s_tet[1]);
7259 nod_pt->set_value(2, s_tet[2]);
7260 tet_el_pt->interpolated_x(s_tet, x_tet);
7261 nod_pt->x(0) = x_tet[0];
7262 nod_pt->x(1) = x_tet[1];
7263 nod_pt->x(2) = x_tet[2];
7266 tet_el_pt->local_coordinate_of_node(6, s_tet);
7267 nod_pt->set_value(0, s_tet[0]);
7268 nod_pt->set_value(1, s_tet[1]);
7269 nod_pt->set_value(2, s_tet[2]);
7270 tet_el_pt->interpolated_x(s_tet, x_tet);
7271 nod_pt->x(0) = x_tet[0];
7272 nod_pt->x(1) = x_tet[1];
7273 nod_pt->x(2) = x_tet[2];
7276 tet_el_pt->local_coordinate_of_node(9, s_tet);
7277 nod_pt->set_value(0, s_tet[0]);
7278 nod_pt->set_value(1, s_tet[1]);
7279 nod_pt->set_value(2, s_tet[2]);
7280 tet_el_pt->interpolated_x(s_tet, x_tet);
7281 nod_pt->x(0) = x_tet[0];
7282 nod_pt->x(1) = x_tet[1];
7283 nod_pt->x(2) = x_tet[2];
7288 s_tet[0] = 1.0 / 3.0;
7289 s_tet[1] = 1.0 / 3.0;
7291 nod_pt->set_value(0, s_tet[0]);
7292 nod_pt->set_value(1, s_tet[1]);
7293 nod_pt->set_value(2, s_tet[2]);
7294 tet_el_pt->interpolated_x(s_tet, x_tet);
7295 nod_pt->x(0) = x_tet[0];
7296 nod_pt->x(1) = x_tet[1];
7297 nod_pt->x(2) = x_tet[2];
7300 tet_el_pt->local_coordinate_of_node(8, s_tet);
7301 nod_pt->set_value(0, s_tet[0]);
7302 nod_pt->set_value(1, s_tet[1]);
7303 nod_pt->set_value(2, s_tet[2]);
7304 tet_el_pt->interpolated_x(s_tet, x_tet);
7305 nod_pt->x(0) = x_tet[0];
7306 nod_pt->x(1) = x_tet[1];
7307 nod_pt->x(2) = x_tet[2];
7312 s_tet[0] = 1.0 / 3.0;
7314 s_tet[2] = 1.0 / 3.0;
7315 nod_pt->set_value(0, s_tet[0]);
7316 nod_pt->set_value(1, s_tet[1]);
7317 nod_pt->set_value(2, s_tet[2]);
7318 tet_el_pt->interpolated_x(s_tet, x_tet);
7319 nod_pt->x(0) = x_tet[0];
7320 nod_pt->x(1) = x_tet[1];
7321 nod_pt->x(2) = x_tet[2];
7327 s_tet[1] = 1.0 / 3.0;
7328 s_tet[2] = 1.0 / 3.0;
7329 nod_pt->set_value(0, s_tet[0]);
7330 nod_pt->set_value(1, s_tet[1]);
7331 nod_pt->set_value(2, s_tet[2]);
7332 tet_el_pt->interpolated_x(s_tet, x_tet);
7333 nod_pt->x(0) = x_tet[0];
7334 nod_pt->x(1) = x_tet[1];
7335 nod_pt->x(2) = x_tet[2];
7342 nod_pt->set_value(0, s_tet[0]);
7343 nod_pt->set_value(1, s_tet[1]);
7344 nod_pt->set_value(2, s_tet[2]);
7345 tet_el_pt->interpolated_x(s_tet, x_tet);
7346 nod_pt->x(0) = x_tet[0];
7347 nod_pt->x(1) = x_tet[1];
7348 nod_pt->x(2) = x_tet[2];
7355 FiniteElement* el_pt =
new ELEMENT;
7356 brick_el2_pt = el_pt;
7360 tet_el_pt->node_pt(3), tet_el_pt->node_pt(0), tet_el_pt->node_pt(2));
7363 tet_el_pt->node_pt(3), tet_el_pt->node_pt(1), tet_el_pt->node_pt(2));
7366 tet_el_pt->node_pt(3), tet_el_pt->node_pt(1), tet_el_pt->node_pt(0));
7369 Vector<Vector<unsigned>> tet_edge_node(3);
7370 tet_edge_node[0].resize(2);
7371 tet_edge_node[0][0] = 6;
7372 tet_edge_node[0][1] = 0;
7373 tet_edge_node[1].resize(2);
7374 tet_edge_node[1][0] = 9;
7375 tet_edge_node[1][1] = 1;
7376 tet_edge_node[2].resize(2);
7377 tet_edge_node[2][0] = 8;
7378 tet_edge_node[2][1] = 2;
7381 unsigned central_tet_vertex = 3;
7383 Node* tet_node_pt = 0;
7384 Node* old_node_pt = 0;
7391 tet_node_pt = tet_el_pt->node_pt(central_tet_vertex);
7392 old_node_pt = tet_node_node_pt[tet_node_pt];
7393 if (old_node_pt == 0)
7395 Node* new_node_pt = 0;
7396 if (tet_node_pt->is_on_boundary())
7398 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7402 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7404 tet_node_node_pt[tet_node_pt] = new_node_pt;
7405 Node_pt.push_back(new_node_pt);
7406 Vector<double>
s(3);
7407 Vector<double> s_tet(3);
7408 Vector<double> x_tet(3);
7409 el_pt->local_coordinate_of_node(
j,
s);
7410 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7411 tet_el_pt->interpolated_x(s_tet, x_tet);
7412 new_node_pt->x(0) = x_tet[0];
7413 new_node_pt->x(1) = x_tet[1];
7414 new_node_pt->x(2) = x_tet[2];
7419 el_pt->node_pt(
j) = old_node_pt;
7429 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[0][0]);
7430 old_node_pt = tet_node_node_pt[tet_node_pt];
7431 if (old_node_pt == 0)
7433 Node* new_node_pt = 0;
7434 if (tet_node_pt->is_on_boundary())
7436 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7440 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7442 tet_node_node_pt[tet_node_pt] = new_node_pt;
7443 Node_pt.push_back(new_node_pt);
7444 Vector<double>
s(3);
7445 Vector<double> s_tet(3);
7446 Vector<double> x_tet(3);
7447 el_pt->local_coordinate_of_node(
j,
s);
7448 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7449 tet_el_pt->interpolated_x(s_tet, x_tet);
7450 new_node_pt->x(0) = x_tet[0];
7451 new_node_pt->x(1) = x_tet[1];
7452 new_node_pt->x(2) = x_tet[2];
7457 el_pt->node_pt(
j) = old_node_pt;
7467 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[1][0]);
7468 old_node_pt = tet_node_node_pt[tet_node_pt];
7469 if (old_node_pt == 0)
7471 Node* new_node_pt = 0;
7472 if (tet_node_pt->is_on_boundary())
7474 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7478 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7480 tet_node_node_pt[tet_node_pt] = new_node_pt;
7481 Node_pt.push_back(new_node_pt);
7482 Vector<double>
s(3);
7483 Vector<double> s_tet(3);
7484 Vector<double> x_tet(3);
7485 el_pt->local_coordinate_of_node(
j,
s);
7486 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7487 tet_el_pt->interpolated_x(s_tet, x_tet);
7488 new_node_pt->x(0) = x_tet[0];
7489 new_node_pt->x(1) = x_tet[1];
7490 new_node_pt->x(2) = x_tet[2];
7495 el_pt->node_pt(
j) = old_node_pt;
7505 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[2][0]);
7506 old_node_pt = tet_node_node_pt[tet_node_pt];
7507 if (old_node_pt == 0)
7509 Node* new_node_pt = 0;
7510 if (tet_node_pt->is_on_boundary())
7512 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7516 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7518 tet_node_node_pt[tet_node_pt] = new_node_pt;
7519 Node_pt.push_back(new_node_pt);
7520 Vector<double>
s(3);
7521 Vector<double> s_tet(3);
7522 Vector<double> x_tet(3);
7523 el_pt->local_coordinate_of_node(
j,
s);
7524 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7525 tet_el_pt->interpolated_x(s_tet, x_tet);
7526 new_node_pt->x(0) = x_tet[0];
7527 new_node_pt->x(1) = x_tet[1];
7528 new_node_pt->x(2) = x_tet[2];
7533 el_pt->node_pt(
j) = old_node_pt;
7543 old_node_pt = tet_face_node_pt[
face0];
7544 if (old_node_pt == 0)
7546 Node* new_node_pt = 0;
7547 if (
face0.is_boundary_face())
7549 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7553 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7555 tet_face_node_pt[
face0] = new_node_pt;
7556 Node_pt.push_back(new_node_pt);
7557 Vector<double>
s(3);
7558 Vector<double> s_tet(3);
7559 Vector<double> x_tet(3);
7560 el_pt->local_coordinate_of_node(
j,
s);
7561 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7562 tet_el_pt->interpolated_x(s_tet, x_tet);
7563 new_node_pt->x(0) = x_tet[0];
7564 new_node_pt->x(1) = x_tet[1];
7565 new_node_pt->x(2) = x_tet[2];
7570 el_pt->node_pt(
j) = old_node_pt;
7579 old_node_pt = tet_face_node_pt[
face1];
7580 if (old_node_pt == 0)
7582 Node* new_node_pt = 0;
7583 if (
face1.is_boundary_face())
7585 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7589 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7591 tet_face_node_pt[
face1] = new_node_pt;
7592 Node_pt.push_back(new_node_pt);
7593 Vector<double>
s(3);
7594 Vector<double> s_tet(3);
7595 Vector<double> x_tet(3);
7596 el_pt->local_coordinate_of_node(
j,
s);
7597 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7598 tet_el_pt->interpolated_x(s_tet, x_tet);
7599 new_node_pt->x(0) = x_tet[0];
7600 new_node_pt->x(1) = x_tet[1];
7601 new_node_pt->x(2) = x_tet[2];
7606 el_pt->node_pt(
j) = old_node_pt;
7615 old_node_pt = tet_face_node_pt[
face2];
7616 if (old_node_pt == 0)
7618 Node* new_node_pt = 0;
7619 if (
face2.is_boundary_face())
7621 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7625 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7627 tet_face_node_pt[
face2] = new_node_pt;
7628 Node_pt.push_back(new_node_pt);
7629 Vector<double>
s(3);
7630 Vector<double> s_tet(3);
7631 Vector<double> x_tet(3);
7632 el_pt->local_coordinate_of_node(
j,
s);
7633 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7634 tet_el_pt->interpolated_x(s_tet, x_tet);
7635 new_node_pt->x(0) = x_tet[0];
7636 new_node_pt->x(1) = x_tet[1];
7637 new_node_pt->x(2) = x_tet[2];
7642 el_pt->node_pt(
j) = old_node_pt;
7653 el_pt->node_pt(
j) = centroid_node_pt;
7663 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7664 Node_pt.push_back(new_node_pt);
7665 Vector<double>
s(3);
7666 Vector<double> s_tet(3);
7667 Vector<double> x_tet(3);
7668 el_pt->local_coordinate_of_node(
j,
s);
7669 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7670 tet_el_pt->interpolated_x(s_tet, x_tet);
7671 new_node_pt->x(0) = x_tet[0];
7672 new_node_pt->x(1) = x_tet[1];
7673 new_node_pt->x(2) = x_tet[2];
7682 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
7683 old_node_pt = brick_edge_node_pt[edge];
7684 if (old_node_pt == 0)
7686 Node* new_node_pt = 0;
7687 if (edge.is_boundary_edge())
7689 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7693 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7695 brick_edge_node_pt[edge] = new_node_pt;
7696 Node_pt.push_back(new_node_pt);
7697 Vector<double>
s(3);
7698 Vector<double> s_tet(3);
7699 Vector<double> x_tet(3);
7700 el_pt->local_coordinate_of_node(
j,
s);
7701 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7702 tet_el_pt->interpolated_x(s_tet, x_tet);
7703 new_node_pt->x(0) = x_tet[0];
7704 new_node_pt->x(1) = x_tet[1];
7705 new_node_pt->x(2) = x_tet[2];
7710 el_pt->node_pt(
j) = old_node_pt;
7720 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
7721 old_node_pt = brick_edge_node_pt[edge];
7722 if (old_node_pt == 0)
7724 Node* new_node_pt = 0;
7725 if (edge.is_boundary_edge())
7727 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7731 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7733 brick_edge_node_pt[edge] = new_node_pt;
7734 Node_pt.push_back(new_node_pt);
7735 Vector<double>
s(3);
7736 Vector<double> s_tet(3);
7737 Vector<double> x_tet(3);
7738 el_pt->local_coordinate_of_node(
j,
s);
7739 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7740 tet_el_pt->interpolated_x(s_tet, x_tet);
7741 new_node_pt->x(0) = x_tet[0];
7742 new_node_pt->x(1) = x_tet[1];
7743 new_node_pt->x(2) = x_tet[2];
7748 el_pt->node_pt(
j) = old_node_pt;
7757 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
7758 old_node_pt = brick_edge_node_pt[edge];
7759 if (old_node_pt == 0)
7761 Node* new_node_pt = 0;
7762 if (edge.is_boundary_edge())
7764 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7768 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7770 brick_edge_node_pt[edge] = new_node_pt;
7771 Node_pt.push_back(new_node_pt);
7772 Vector<double>
s(3);
7773 Vector<double> s_tet(3);
7774 Vector<double> x_tet(3);
7775 el_pt->local_coordinate_of_node(
j,
s);
7776 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7777 tet_el_pt->interpolated_x(s_tet, x_tet);
7778 new_node_pt->x(0) = x_tet[0];
7779 new_node_pt->x(1) = x_tet[1];
7780 new_node_pt->x(2) = x_tet[2];
7785 el_pt->node_pt(
j) = old_node_pt;
7794 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
7795 old_node_pt = brick_edge_node_pt[edge];
7796 if (old_node_pt == 0)
7798 Node* new_node_pt = 0;
7799 if (edge.is_boundary_edge())
7801 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7805 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7807 brick_edge_node_pt[edge] = new_node_pt;
7808 Node_pt.push_back(new_node_pt);
7809 Vector<double>
s(3);
7810 Vector<double> s_tet(3);
7811 Vector<double> x_tet(3);
7812 el_pt->local_coordinate_of_node(
j,
s);
7813 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7814 tet_el_pt->interpolated_x(s_tet, x_tet);
7815 new_node_pt->x(0) = x_tet[0];
7816 new_node_pt->x(1) = x_tet[1];
7817 new_node_pt->x(2) = x_tet[2];
7822 el_pt->node_pt(
j) = old_node_pt;
7831 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
7832 old_node_pt = brick_edge_node_pt[edge];
7833 if (old_node_pt == 0)
7835 Node* new_node_pt = 0;
7836 if (edge.is_boundary_edge())
7838 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7842 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7844 brick_edge_node_pt[edge] = new_node_pt;
7845 Node_pt.push_back(new_node_pt);
7846 Vector<double>
s(3);
7847 Vector<double> s_tet(3);
7848 Vector<double> x_tet(3);
7849 el_pt->local_coordinate_of_node(
j,
s);
7850 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7851 tet_el_pt->interpolated_x(s_tet, x_tet);
7852 new_node_pt->x(0) = x_tet[0];
7853 new_node_pt->x(1) = x_tet[1];
7854 new_node_pt->x(2) = x_tet[2];
7859 el_pt->node_pt(
j) = old_node_pt;
7869 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
7870 old_node_pt = brick_edge_node_pt[edge];
7871 if (old_node_pt == 0)
7873 Node* new_node_pt = 0;
7874 if (edge.is_boundary_edge())
7876 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7880 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7882 brick_edge_node_pt[edge] = new_node_pt;
7883 Node_pt.push_back(new_node_pt);
7884 Vector<double>
s(3);
7885 Vector<double> s_tet(3);
7886 Vector<double> x_tet(3);
7887 el_pt->local_coordinate_of_node(
j,
s);
7888 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7889 tet_el_pt->interpolated_x(s_tet, x_tet);
7890 new_node_pt->x(0) = x_tet[0];
7891 new_node_pt->x(1) = x_tet[1];
7892 new_node_pt->x(2) = x_tet[2];
7897 el_pt->node_pt(
j) = old_node_pt;
7906 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
7907 old_node_pt = brick_edge_node_pt[edge];
7908 if (old_node_pt == 0)
7910 Node* new_node_pt = 0;
7911 if (edge.is_boundary_edge())
7913 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7917 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7919 brick_edge_node_pt[edge] = new_node_pt;
7920 Node_pt.push_back(new_node_pt);
7921 Vector<double>
s(3);
7922 Vector<double> s_tet(3);
7923 Vector<double> x_tet(3);
7924 el_pt->local_coordinate_of_node(
j,
s);
7925 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7926 tet_el_pt->interpolated_x(s_tet, x_tet);
7927 new_node_pt->x(0) = x_tet[0];
7928 new_node_pt->x(1) = x_tet[1];
7929 new_node_pt->x(2) = x_tet[2];
7934 el_pt->node_pt(
j) = old_node_pt;
7944 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
7945 old_node_pt = brick_edge_node_pt[edge];
7946 if (old_node_pt == 0)
7948 Node* new_node_pt = 0;
7949 if (edge.is_boundary_edge())
7951 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7955 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7957 brick_edge_node_pt[edge] = new_node_pt;
7958 Node_pt.push_back(new_node_pt);
7959 Vector<double>
s(3);
7960 Vector<double> s_tet(3);
7961 Vector<double> x_tet(3);
7962 el_pt->local_coordinate_of_node(
j,
s);
7963 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
7964 tet_el_pt->interpolated_x(s_tet, x_tet);
7965 new_node_pt->x(0) = x_tet[0];
7966 new_node_pt->x(1) = x_tet[1];
7967 new_node_pt->x(2) = x_tet[2];
7972 el_pt->node_pt(
j) = old_node_pt;
7981 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
7982 old_node_pt = brick_edge_node_pt[edge];
7983 if (old_node_pt == 0)
7985 Node* new_node_pt = 0;
7986 if (edge.is_boundary_edge())
7988 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
7992 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
7994 brick_edge_node_pt[edge] = new_node_pt;
7995 Node_pt.push_back(new_node_pt);
7996 Vector<double>
s(3);
7997 Vector<double> s_tet(3);
7998 Vector<double> x_tet(3);
7999 el_pt->local_coordinate_of_node(
j,
s);
8000 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
8001 tet_el_pt->interpolated_x(s_tet, x_tet);
8002 new_node_pt->x(0) = x_tet[0];
8003 new_node_pt->x(1) = x_tet[1];
8004 new_node_pt->x(2) = x_tet[2];
8009 el_pt->node_pt(
j) = old_node_pt;
8019 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
8020 old_node_pt = brick_edge_node_pt[edge];
8021 if (old_node_pt == 0)
8023 Node* new_node_pt = 0;
8024 if (edge.is_boundary_edge())
8026 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8030 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8032 brick_edge_node_pt[edge] = new_node_pt;
8033 Node_pt.push_back(new_node_pt);
8034 Vector<double>
s(3);
8035 Vector<double> s_tet(3);
8036 Vector<double> x_tet(3);
8037 el_pt->local_coordinate_of_node(
j,
s);
8038 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
8039 tet_el_pt->interpolated_x(s_tet, x_tet);
8040 new_node_pt->x(0) = x_tet[0];
8041 new_node_pt->x(1) = x_tet[1];
8042 new_node_pt->x(2) = x_tet[2];
8047 el_pt->node_pt(
j) = old_node_pt;
8057 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
8058 old_node_pt = brick_edge_node_pt[edge];
8059 if (old_node_pt == 0)
8061 Node* new_node_pt = 0;
8062 if (edge.is_boundary_edge())
8064 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8068 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8070 brick_edge_node_pt[edge] = new_node_pt;
8071 Node_pt.push_back(new_node_pt);
8072 Vector<double>
s(3);
8073 Vector<double> s_tet(3);
8074 Vector<double> x_tet(3);
8075 el_pt->local_coordinate_of_node(
j,
s);
8076 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
8077 tet_el_pt->interpolated_x(s_tet, x_tet);
8078 new_node_pt->x(0) = x_tet[0];
8079 new_node_pt->x(1) = x_tet[1];
8080 new_node_pt->x(2) = x_tet[2];
8085 el_pt->node_pt(
j) = old_node_pt;
8095 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
8096 old_node_pt = brick_edge_node_pt[edge];
8097 if (old_node_pt == 0)
8099 Node* new_node_pt = 0;
8100 if (edge.is_boundary_edge())
8102 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8106 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8108 brick_edge_node_pt[edge] = new_node_pt;
8109 Node_pt.push_back(new_node_pt);
8110 Vector<double>
s(3);
8111 Vector<double> s_tet(3);
8112 Vector<double> x_tet(3);
8113 el_pt->local_coordinate_of_node(
j,
s);
8114 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
8115 tet_el_pt->interpolated_x(s_tet, x_tet);
8116 new_node_pt->x(0) = x_tet[0];
8117 new_node_pt->x(1) = x_tet[1];
8118 new_node_pt->x(2) = x_tet[2];
8123 el_pt->node_pt(
j) = old_node_pt;
8134 TFace face(tet_el_pt->node_pt(central_tet_vertex),
8135 tet_el_pt->node_pt(tet_edge_node[0][0]),
8136 tet_el_pt->node_pt(tet_edge_node[2][0]));
8138 old_node_pt = tet_face_node_pt[face];
8139 if (old_node_pt == 0)
8141 Node* new_node_pt = 0;
8142 if (face.is_boundary_face())
8144 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8148 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8150 tet_face_node_pt[face] = new_node_pt;
8151 Node_pt.push_back(new_node_pt);
8152 Vector<double>
s(3);
8153 Vector<double> s_tet(3);
8154 Vector<double> x_tet(3);
8155 el_pt->local_coordinate_of_node(
j,
s);
8156 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
8157 tet_el_pt->interpolated_x(s_tet, x_tet);
8158 new_node_pt->x(0) = x_tet[0];
8159 new_node_pt->x(1) = x_tet[1];
8160 new_node_pt->x(2) = x_tet[2];
8165 el_pt->node_pt(
j) = old_node_pt;
8176 TFace face(tet_el_pt->node_pt(central_tet_vertex),
8177 tet_el_pt->node_pt(tet_edge_node[1][0]),
8178 tet_el_pt->node_pt(tet_edge_node[2][0]));
8179 old_node_pt = tet_face_node_pt[face];
8180 if (old_node_pt == 0)
8182 Node* new_node_pt = 0;
8183 if (face.is_boundary_face())
8185 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8189 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8191 tet_face_node_pt[face] = new_node_pt;
8192 Node_pt.push_back(new_node_pt);
8193 Vector<double>
s(3);
8194 Vector<double> s_tet(3);
8195 Vector<double> x_tet(3);
8196 el_pt->local_coordinate_of_node(
j,
s);
8197 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
8198 tet_el_pt->interpolated_x(s_tet, x_tet);
8199 new_node_pt->x(0) = x_tet[0];
8200 new_node_pt->x(1) = x_tet[1];
8201 new_node_pt->x(2) = x_tet[2];
8206 el_pt->node_pt(
j) = old_node_pt;
8217 TFace face(tet_el_pt->node_pt(central_tet_vertex),
8218 tet_el_pt->node_pt(tet_edge_node[0][0]),
8219 tet_el_pt->node_pt(tet_edge_node[1][0]));
8221 old_node_pt = tet_face_node_pt[face];
8222 if (old_node_pt == 0)
8224 Node* new_node_pt = 0;
8225 if (face.is_boundary_face())
8227 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8231 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8233 tet_face_node_pt[face] = new_node_pt;
8234 Node_pt.push_back(new_node_pt);
8235 Vector<double>
s(3);
8236 Vector<double> s_tet(3);
8237 Vector<double> x_tet(3);
8238 el_pt->local_coordinate_of_node(
j,
s);
8239 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
8240 tet_el_pt->interpolated_x(s_tet, x_tet);
8241 new_node_pt->x(0) = x_tet[0];
8242 new_node_pt->x(1) = x_tet[1];
8243 new_node_pt->x(2) = x_tet[2];
8248 el_pt->node_pt(
j) = old_node_pt;
8256 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8257 Node_pt.push_back(new_node_pt);
8258 Vector<double>
s(3);
8259 Vector<double> s_tet(3);
8260 Vector<double> x_tet(3);
8261 el_pt->local_coordinate_of_node(
j,
s);
8262 dummy_q_el_pt[2]->interpolated_s_tet(
s, s_tet);
8263 top_mid_face_node2_pt = new_node_pt;
8264 tet_el_pt->interpolated_x(s_tet, x_tet);
8265 new_node_pt->x(0) = x_tet[0];
8266 new_node_pt->x(1) = x_tet[1];
8267 new_node_pt->x(2) = x_tet[2];
8276 el_pt->node_pt(
j) = back_mid_face_node0_pt;
8285 el_pt->node_pt(
j) = right_mid_face_node1_pt;
8294 for (
unsigned j = 0;
j < 8;
j++)
8296 Node* nod_pt = dummy_q_el_pt[3]->node_pt(
j);
8297 Vector<double> s_tet(3);
8298 Vector<double> x_tet(3);
8302 tet_el_pt->local_coordinate_of_node(2, s_tet);
8303 nod_pt->set_value(0, s_tet[0]);
8304 nod_pt->set_value(1, s_tet[1]);
8305 nod_pt->set_value(2, s_tet[2]);
8306 tet_el_pt->interpolated_x(s_tet, x_tet);
8307 nod_pt->x(0) = x_tet[0];
8308 nod_pt->x(1) = x_tet[1];
8309 nod_pt->x(2) = x_tet[2];
8312 tet_el_pt->local_coordinate_of_node(7, s_tet);
8313 nod_pt->set_value(0, s_tet[0]);
8314 nod_pt->set_value(1, s_tet[1]);
8315 nod_pt->set_value(2, s_tet[2]);
8316 tet_el_pt->interpolated_x(s_tet, x_tet);
8317 nod_pt->x(0) = x_tet[0];
8318 nod_pt->x(1) = x_tet[1];
8319 nod_pt->x(2) = x_tet[2];
8322 tet_el_pt->local_coordinate_of_node(5, s_tet);
8323 nod_pt->set_value(0, s_tet[0]);
8324 nod_pt->set_value(1, s_tet[1]);
8325 nod_pt->set_value(2, s_tet[2]);
8326 tet_el_pt->interpolated_x(s_tet, x_tet);
8327 nod_pt->x(0) = x_tet[0];
8328 nod_pt->x(1) = x_tet[1];
8329 nod_pt->x(2) = x_tet[2];
8334 s_tet[0] = 1.0 / 3.0;
8335 s_tet[1] = 1.0 / 3.0;
8336 s_tet[2] = 1.0 / 3.0;
8337 nod_pt->set_value(0, s_tet[0]);
8338 nod_pt->set_value(1, s_tet[1]);
8339 nod_pt->set_value(2, s_tet[2]);
8340 tet_el_pt->interpolated_x(s_tet, x_tet);
8341 nod_pt->x(0) = x_tet[0];
8342 nod_pt->x(1) = x_tet[1];
8343 nod_pt->x(2) = x_tet[2];
8346 tet_el_pt->local_coordinate_of_node(8, s_tet);
8347 nod_pt->set_value(0, s_tet[0]);
8348 nod_pt->set_value(1, s_tet[1]);
8349 nod_pt->set_value(2, s_tet[2]);
8350 tet_el_pt->interpolated_x(s_tet, x_tet);
8351 nod_pt->x(0) = x_tet[0];
8352 nod_pt->x(1) = x_tet[1];
8353 nod_pt->x(2) = x_tet[2];
8359 s_tet[1] = 1.0 / 3.0;
8360 s_tet[2] = 1.0 / 3.0;
8361 nod_pt->set_value(0, s_tet[0]);
8362 nod_pt->set_value(1, s_tet[1]);
8363 nod_pt->set_value(2, s_tet[2]);
8364 tet_el_pt->interpolated_x(s_tet, x_tet);
8365 nod_pt->x(0) = x_tet[0];
8366 nod_pt->x(1) = x_tet[1];
8367 nod_pt->x(2) = x_tet[2];
8372 s_tet[0] = 1.0 / 3.0;
8374 s_tet[2] = 1.0 / 3.0;
8375 nod_pt->set_value(0, s_tet[0]);
8376 nod_pt->set_value(1, s_tet[1]);
8377 nod_pt->set_value(2, s_tet[2]);
8378 tet_el_pt->interpolated_x(s_tet, x_tet);
8379 nod_pt->x(0) = x_tet[0];
8380 nod_pt->x(1) = x_tet[1];
8381 nod_pt->x(2) = x_tet[2];
8388 nod_pt->set_value(0, s_tet[0]);
8389 nod_pt->set_value(1, s_tet[1]);
8390 nod_pt->set_value(2, s_tet[2]);
8391 tet_el_pt->interpolated_x(s_tet, x_tet);
8392 nod_pt->x(0) = x_tet[0];
8393 nod_pt->x(1) = x_tet[1];
8394 nod_pt->x(2) = x_tet[2];
8401 FiniteElement* el_pt =
new ELEMENT;
8402 brick_el3_pt = el_pt;
8406 tet_el_pt->node_pt(1), tet_el_pt->node_pt(2), tet_el_pt->node_pt(3));
8409 tet_el_pt->node_pt(0), tet_el_pt->node_pt(2), tet_el_pt->node_pt(3));
8412 tet_el_pt->node_pt(0), tet_el_pt->node_pt(1), tet_el_pt->node_pt(2));
8415 Vector<Vector<unsigned>> tet_edge_node(3);
8416 tet_edge_node[0].resize(2);
8417 tet_edge_node[0][0] = 7;
8418 tet_edge_node[0][1] = 1;
8419 tet_edge_node[1].resize(2);
8420 tet_edge_node[1][0] = 5;
8421 tet_edge_node[1][1] = 0;
8422 tet_edge_node[2].resize(2);
8423 tet_edge_node[2][0] = 8;
8424 tet_edge_node[2][1] = 3;
8427 unsigned central_tet_vertex = 2;
8429 Node* tet_node_pt = 0;
8430 Node* old_node_pt = 0;
8437 tet_node_pt = tet_el_pt->node_pt(central_tet_vertex);
8438 old_node_pt = tet_node_node_pt[tet_node_pt];
8439 if (old_node_pt == 0)
8441 Node* new_node_pt = 0;
8442 if (tet_node_pt->is_on_boundary())
8444 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8448 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8450 tet_node_node_pt[tet_node_pt] = new_node_pt;
8451 Node_pt.push_back(new_node_pt);
8452 Vector<double>
s(3);
8453 Vector<double> s_tet(3);
8454 Vector<double> x_tet(3);
8455 el_pt->local_coordinate_of_node(
j,
s);
8456 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8457 tet_el_pt->interpolated_x(s_tet, x_tet);
8458 new_node_pt->x(0) = x_tet[0];
8459 new_node_pt->x(1) = x_tet[1];
8460 new_node_pt->x(2) = x_tet[2];
8465 el_pt->node_pt(
j) = old_node_pt;
8475 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[0][0]);
8476 old_node_pt = tet_node_node_pt[tet_node_pt];
8477 if (old_node_pt == 0)
8479 Node* new_node_pt = 0;
8480 if (tet_node_pt->is_on_boundary())
8482 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8486 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8488 tet_node_node_pt[tet_node_pt] = new_node_pt;
8489 Node_pt.push_back(new_node_pt);
8490 Vector<double>
s(3);
8491 Vector<double> s_tet(3);
8492 Vector<double> x_tet(3);
8493 el_pt->local_coordinate_of_node(
j,
s);
8494 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8495 tet_el_pt->interpolated_x(s_tet, x_tet);
8496 new_node_pt->x(0) = x_tet[0];
8497 new_node_pt->x(1) = x_tet[1];
8498 new_node_pt->x(2) = x_tet[2];
8503 el_pt->node_pt(
j) = old_node_pt;
8513 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[1][0]);
8514 old_node_pt = tet_node_node_pt[tet_node_pt];
8515 if (old_node_pt == 0)
8517 Node* new_node_pt = 0;
8518 if (tet_node_pt->is_on_boundary())
8520 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8524 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8526 tet_node_node_pt[tet_node_pt] = new_node_pt;
8527 Node_pt.push_back(new_node_pt);
8528 Vector<double>
s(3);
8529 Vector<double> s_tet(3);
8530 Vector<double> x_tet(3);
8531 el_pt->local_coordinate_of_node(
j,
s);
8532 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8533 tet_el_pt->interpolated_x(s_tet, x_tet);
8534 new_node_pt->x(0) = x_tet[0];
8535 new_node_pt->x(1) = x_tet[1];
8536 new_node_pt->x(2) = x_tet[2];
8541 el_pt->node_pt(
j) = old_node_pt;
8551 tet_node_pt = tet_el_pt->node_pt(tet_edge_node[2][0]);
8552 old_node_pt = tet_node_node_pt[tet_node_pt];
8553 if (old_node_pt == 0)
8555 Node* new_node_pt = 0;
8556 if (tet_node_pt->is_on_boundary())
8558 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8562 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8564 tet_node_node_pt[tet_node_pt] = new_node_pt;
8565 Node_pt.push_back(new_node_pt);
8566 Vector<double>
s(3);
8567 Vector<double> s_tet(3);
8568 Vector<double> x_tet(3);
8569 el_pt->local_coordinate_of_node(
j,
s);
8570 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8571 tet_el_pt->interpolated_x(s_tet, x_tet);
8572 new_node_pt->x(0) = x_tet[0];
8573 new_node_pt->x(1) = x_tet[1];
8574 new_node_pt->x(2) = x_tet[2];
8579 el_pt->node_pt(
j) = old_node_pt;
8589 old_node_pt = tet_face_node_pt[
face0];
8590 if (old_node_pt == 0)
8592 Node* new_node_pt = 0;
8593 if (
face0.is_boundary_face())
8595 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8599 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8601 tet_face_node_pt[
face0] = new_node_pt;
8602 Node_pt.push_back(new_node_pt);
8603 Vector<double>
s(3);
8604 Vector<double> s_tet(3);
8605 Vector<double> x_tet(3);
8606 el_pt->local_coordinate_of_node(
j,
s);
8607 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8608 tet_el_pt->interpolated_x(s_tet, x_tet);
8609 new_node_pt->x(0) = x_tet[0];
8610 new_node_pt->x(1) = x_tet[1];
8611 new_node_pt->x(2) = x_tet[2];
8616 el_pt->node_pt(
j) = old_node_pt;
8625 old_node_pt = tet_face_node_pt[
face1];
8626 if (old_node_pt == 0)
8628 Node* new_node_pt = 0;
8629 if (
face1.is_boundary_face())
8631 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8635 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8637 tet_face_node_pt[
face1] = new_node_pt;
8638 Node_pt.push_back(new_node_pt);
8639 Vector<double>
s(3);
8640 Vector<double> s_tet(3);
8641 Vector<double> x_tet(3);
8642 el_pt->local_coordinate_of_node(
j,
s);
8643 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8644 tet_el_pt->interpolated_x(s_tet, x_tet);
8645 new_node_pt->x(0) = x_tet[0];
8646 new_node_pt->x(1) = x_tet[1];
8647 new_node_pt->x(2) = x_tet[2];
8652 el_pt->node_pt(
j) = old_node_pt;
8661 old_node_pt = tet_face_node_pt[
face2];
8662 if (old_node_pt == 0)
8664 Node* new_node_pt = 0;
8665 if (
face2.is_boundary_face())
8667 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8671 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8673 tet_face_node_pt[
face2] = new_node_pt;
8674 Node_pt.push_back(new_node_pt);
8675 Vector<double>
s(3);
8676 Vector<double> s_tet(3);
8677 Vector<double> x_tet(3);
8678 el_pt->local_coordinate_of_node(
j,
s);
8679 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8680 tet_el_pt->interpolated_x(s_tet, x_tet);
8681 new_node_pt->x(0) = x_tet[0];
8682 new_node_pt->x(1) = x_tet[1];
8683 new_node_pt->x(2) = x_tet[2];
8688 el_pt->node_pt(
j) = old_node_pt;
8699 el_pt->node_pt(
j) = centroid_node_pt;
8709 Node* new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8710 Node_pt.push_back(new_node_pt);
8711 Vector<double>
s(3);
8712 Vector<double> s_tet(3);
8713 Vector<double> x_tet(3);
8714 el_pt->local_coordinate_of_node(
j,
s);
8715 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8716 tet_el_pt->interpolated_x(s_tet, x_tet);
8717 new_node_pt->x(0) = x_tet[0];
8718 new_node_pt->x(1) = x_tet[1];
8719 new_node_pt->x(2) = x_tet[2];
8728 Edge edge(el_pt->node_pt(0), el_pt->node_pt(2));
8729 old_node_pt = brick_edge_node_pt[edge];
8730 if (old_node_pt == 0)
8732 Node* new_node_pt = 0;
8733 if (edge.is_boundary_edge())
8735 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8739 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8741 brick_edge_node_pt[edge] = new_node_pt;
8742 Node_pt.push_back(new_node_pt);
8743 Vector<double>
s(3);
8744 Vector<double> s_tet(3);
8745 Vector<double> x_tet(3);
8746 el_pt->local_coordinate_of_node(
j,
s);
8747 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8748 tet_el_pt->interpolated_x(s_tet, x_tet);
8749 new_node_pt->x(0) = x_tet[0];
8750 new_node_pt->x(1) = x_tet[1];
8751 new_node_pt->x(2) = x_tet[2];
8756 el_pt->node_pt(
j) = old_node_pt;
8766 Edge edge(el_pt->node_pt(0), el_pt->node_pt(6));
8767 old_node_pt = brick_edge_node_pt[edge];
8768 if (old_node_pt == 0)
8770 Node* new_node_pt = 0;
8771 if (edge.is_boundary_edge())
8773 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8777 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8779 brick_edge_node_pt[edge] = new_node_pt;
8780 Node_pt.push_back(new_node_pt);
8781 Vector<double>
s(3);
8782 Vector<double> s_tet(3);
8783 Vector<double> x_tet(3);
8784 el_pt->local_coordinate_of_node(
j,
s);
8785 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8786 tet_el_pt->interpolated_x(s_tet, x_tet);
8787 new_node_pt->x(0) = x_tet[0];
8788 new_node_pt->x(1) = x_tet[1];
8789 new_node_pt->x(2) = x_tet[2];
8794 el_pt->node_pt(
j) = old_node_pt;
8803 Edge edge(el_pt->node_pt(2), el_pt->node_pt(8));
8804 old_node_pt = brick_edge_node_pt[edge];
8805 if (old_node_pt == 0)
8807 Node* new_node_pt = 0;
8808 if (edge.is_boundary_edge())
8810 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8814 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8816 brick_edge_node_pt[edge] = new_node_pt;
8817 Node_pt.push_back(new_node_pt);
8818 Vector<double>
s(3);
8819 Vector<double> s_tet(3);
8820 Vector<double> x_tet(3);
8821 el_pt->local_coordinate_of_node(
j,
s);
8822 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8823 tet_el_pt->interpolated_x(s_tet, x_tet);
8824 new_node_pt->x(0) = x_tet[0];
8825 new_node_pt->x(1) = x_tet[1];
8826 new_node_pt->x(2) = x_tet[2];
8831 el_pt->node_pt(
j) = old_node_pt;
8840 Edge edge(el_pt->node_pt(6), el_pt->node_pt(8));
8841 old_node_pt = brick_edge_node_pt[edge];
8842 if (old_node_pt == 0)
8844 Node* new_node_pt = 0;
8845 if (edge.is_boundary_edge())
8847 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8851 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8853 brick_edge_node_pt[edge] = new_node_pt;
8854 Node_pt.push_back(new_node_pt);
8855 Vector<double>
s(3);
8856 Vector<double> s_tet(3);
8857 Vector<double> x_tet(3);
8858 el_pt->local_coordinate_of_node(
j,
s);
8859 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8860 tet_el_pt->interpolated_x(s_tet, x_tet);
8861 new_node_pt->x(0) = x_tet[0];
8862 new_node_pt->x(1) = x_tet[1];
8863 new_node_pt->x(2) = x_tet[2];
8868 el_pt->node_pt(
j) = old_node_pt;
8877 Edge edge(el_pt->node_pt(18), el_pt->node_pt(20));
8878 old_node_pt = brick_edge_node_pt[edge];
8879 if (old_node_pt == 0)
8881 Node* new_node_pt = 0;
8882 if (edge.is_boundary_edge())
8884 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8888 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8890 brick_edge_node_pt[edge] = new_node_pt;
8891 Node_pt.push_back(new_node_pt);
8892 Vector<double>
s(3);
8893 Vector<double> s_tet(3);
8894 Vector<double> x_tet(3);
8895 el_pt->local_coordinate_of_node(
j,
s);
8896 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8897 tet_el_pt->interpolated_x(s_tet, x_tet);
8898 new_node_pt->x(0) = x_tet[0];
8899 new_node_pt->x(1) = x_tet[1];
8900 new_node_pt->x(2) = x_tet[2];
8905 el_pt->node_pt(
j) = old_node_pt;
8915 Edge edge(el_pt->node_pt(18), el_pt->node_pt(24));
8916 old_node_pt = brick_edge_node_pt[edge];
8917 if (old_node_pt == 0)
8919 Node* new_node_pt = 0;
8920 if (edge.is_boundary_edge())
8922 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8926 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8928 brick_edge_node_pt[edge] = new_node_pt;
8929 Node_pt.push_back(new_node_pt);
8930 Vector<double>
s(3);
8931 Vector<double> s_tet(3);
8932 Vector<double> x_tet(3);
8933 el_pt->local_coordinate_of_node(
j,
s);
8934 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8935 tet_el_pt->interpolated_x(s_tet, x_tet);
8936 new_node_pt->x(0) = x_tet[0];
8937 new_node_pt->x(1) = x_tet[1];
8938 new_node_pt->x(2) = x_tet[2];
8943 el_pt->node_pt(
j) = old_node_pt;
8952 Edge edge(el_pt->node_pt(20), el_pt->node_pt(26));
8953 old_node_pt = brick_edge_node_pt[edge];
8954 if (old_node_pt == 0)
8956 Node* new_node_pt = 0;
8957 if (edge.is_boundary_edge())
8959 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
8963 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
8965 brick_edge_node_pt[edge] = new_node_pt;
8966 Node_pt.push_back(new_node_pt);
8967 Vector<double>
s(3);
8968 Vector<double> s_tet(3);
8969 Vector<double> x_tet(3);
8970 el_pt->local_coordinate_of_node(
j,
s);
8971 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
8972 tet_el_pt->interpolated_x(s_tet, x_tet);
8973 new_node_pt->x(0) = x_tet[0];
8974 new_node_pt->x(1) = x_tet[1];
8975 new_node_pt->x(2) = x_tet[2];
8980 el_pt->node_pt(
j) = old_node_pt;
8990 Edge edge(el_pt->node_pt(24), el_pt->node_pt(26));
8991 old_node_pt = brick_edge_node_pt[edge];
8992 if (old_node_pt == 0)
8994 Node* new_node_pt = 0;
8995 if (edge.is_boundary_edge())
8997 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
9001 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
9003 brick_edge_node_pt[edge] = new_node_pt;
9004 Node_pt.push_back(new_node_pt);
9005 Vector<double>
s(3);
9006 Vector<double> s_tet(3);
9007 Vector<double> x_tet(3);
9008 el_pt->local_coordinate_of_node(
j,
s);
9009 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
9010 tet_el_pt->interpolated_x(s_tet, x_tet);
9011 new_node_pt->x(0) = x_tet[0];
9012 new_node_pt->x(1) = x_tet[1];
9013 new_node_pt->x(2) = x_tet[2];
9018 el_pt->node_pt(
j) = old_node_pt;
9027 Edge edge(el_pt->node_pt(0), el_pt->node_pt(18));
9028 old_node_pt = brick_edge_node_pt[edge];
9029 if (old_node_pt == 0)
9031 Node* new_node_pt = 0;
9032 if (edge.is_boundary_edge())
9034 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
9038 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
9040 brick_edge_node_pt[edge] = new_node_pt;
9041 Node_pt.push_back(new_node_pt);
9042 Vector<double>
s(3);
9043 Vector<double> s_tet(3);
9044 Vector<double> x_tet(3);
9045 el_pt->local_coordinate_of_node(
j,
s);
9046 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
9047 tet_el_pt->interpolated_x(s_tet, x_tet);
9048 new_node_pt->x(0) = x_tet[0];
9049 new_node_pt->x(1) = x_tet[1];
9050 new_node_pt->x(2) = x_tet[2];
9055 el_pt->node_pt(
j) = old_node_pt;
9065 Edge edge(el_pt->node_pt(2), el_pt->node_pt(20));
9066 old_node_pt = brick_edge_node_pt[edge];
9067 if (old_node_pt == 0)
9069 Node* new_node_pt = 0;
9070 if (edge.is_boundary_edge())
9072 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
9076 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
9078 brick_edge_node_pt[edge] = new_node_pt;
9079 Node_pt.push_back(new_node_pt);
9080 Vector<double>
s(3);
9081 Vector<double> s_tet(3);
9082 Vector<double> x_tet(3);
9083 el_pt->local_coordinate_of_node(
j,
s);
9084 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
9085 tet_el_pt->interpolated_x(s_tet, x_tet);
9086 new_node_pt->x(0) = x_tet[0];
9087 new_node_pt->x(1) = x_tet[1];
9088 new_node_pt->x(2) = x_tet[2];
9093 el_pt->node_pt(
j) = old_node_pt;
9103 Edge edge(el_pt->node_pt(6), el_pt->node_pt(24));
9104 old_node_pt = brick_edge_node_pt[edge];
9105 if (old_node_pt == 0)
9107 Node* new_node_pt = 0;
9108 if (edge.is_boundary_edge())
9110 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
9114 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
9116 brick_edge_node_pt[edge] = new_node_pt;
9117 Node_pt.push_back(new_node_pt);
9118 Vector<double>
s(3);
9119 Vector<double> s_tet(3);
9120 Vector<double> x_tet(3);
9121 el_pt->local_coordinate_of_node(
j,
s);
9122 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
9123 tet_el_pt->interpolated_x(s_tet, x_tet);
9124 new_node_pt->x(0) = x_tet[0];
9125 new_node_pt->x(1) = x_tet[1];
9126 new_node_pt->x(2) = x_tet[2];
9131 el_pt->node_pt(
j) = old_node_pt;
9141 Edge edge(el_pt->node_pt(8), el_pt->node_pt(26));
9142 old_node_pt = brick_edge_node_pt[edge];
9143 if (old_node_pt == 0)
9145 Node* new_node_pt = 0;
9146 if (edge.is_boundary_edge())
9148 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
9152 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
9154 brick_edge_node_pt[edge] = new_node_pt;
9155 Node_pt.push_back(new_node_pt);
9156 Vector<double>
s(3);
9157 Vector<double> s_tet(3);
9158 Vector<double> x_tet(3);
9159 el_pt->local_coordinate_of_node(
j,
s);
9160 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
9161 tet_el_pt->interpolated_x(s_tet, x_tet);
9162 new_node_pt->x(0) = x_tet[0];
9163 new_node_pt->x(1) = x_tet[1];
9164 new_node_pt->x(2) = x_tet[2];
9169 el_pt->node_pt(
j) = old_node_pt;
9180 TFace face(tet_el_pt->node_pt(central_tet_vertex),
9181 tet_el_pt->node_pt(tet_edge_node[0][0]),
9182 tet_el_pt->node_pt(tet_edge_node[2][0]));
9184 old_node_pt = tet_face_node_pt[face];
9185 if (old_node_pt == 0)
9187 Node* new_node_pt = 0;
9188 if (face.is_boundary_face())
9190 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
9194 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
9196 tet_face_node_pt[face] = new_node_pt;
9197 Node_pt.push_back(new_node_pt);
9198 Vector<double>
s(3);
9199 Vector<double> s_tet(3);
9200 Vector<double> x_tet(3);
9201 el_pt->local_coordinate_of_node(
j,
s);
9202 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
9203 tet_el_pt->interpolated_x(s_tet, x_tet);
9204 new_node_pt->x(0) = x_tet[0];
9205 new_node_pt->x(1) = x_tet[1];
9206 new_node_pt->x(2) = x_tet[2];
9211 el_pt->node_pt(
j) = old_node_pt;
9222 TFace face(tet_el_pt->node_pt(central_tet_vertex),
9223 tet_el_pt->node_pt(tet_edge_node[1][0]),
9224 tet_el_pt->node_pt(tet_edge_node[2][0]));
9226 old_node_pt = tet_face_node_pt[face];
9227 if (old_node_pt == 0)
9229 Node* new_node_pt = 0;
9230 if (face.is_boundary_face())
9232 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
9236 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
9238 tet_face_node_pt[face] = new_node_pt;
9239 Node_pt.push_back(new_node_pt);
9240 Vector<double>
s(3);
9241 Vector<double> s_tet(3);
9242 Vector<double> x_tet(3);
9243 el_pt->local_coordinate_of_node(
j,
s);
9244 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
9245 tet_el_pt->interpolated_x(s_tet, x_tet);
9246 new_node_pt->x(0) = x_tet[0];
9247 new_node_pt->x(1) = x_tet[1];
9248 new_node_pt->x(2) = x_tet[2];
9253 el_pt->node_pt(
j) = old_node_pt;
9264 TFace face(tet_el_pt->node_pt(central_tet_vertex),
9265 tet_el_pt->node_pt(tet_edge_node[0][0]),
9266 tet_el_pt->node_pt(tet_edge_node[1][0]));
9268 old_node_pt = tet_face_node_pt[face];
9269 if (old_node_pt == 0)
9271 Node* new_node_pt = 0;
9272 if (face.is_boundary_face())
9274 new_node_pt = el_pt->construct_boundary_node(
j, time_stepper_pt);
9278 new_node_pt = el_pt->construct_node(
j, time_stepper_pt);
9280 tet_face_node_pt[face] = new_node_pt;
9281 Node_pt.push_back(new_node_pt);
9282 Vector<double>
s(3);
9283 Vector<double> s_tet(3);
9284 Vector<double> x_tet(3);
9285 el_pt->local_coordinate_of_node(
j,
s);
9286 dummy_q_el_pt[3]->interpolated_s_tet(
s, s_tet);
9287 tet_el_pt->interpolated_x(s_tet, x_tet);
9288 new_node_pt->x(0) = x_tet[0];
9289 new_node_pt->x(1) = x_tet[1];
9290 new_node_pt->x(2) = x_tet[2];
9295 el_pt->node_pt(
j) = old_node_pt;
9305 el_pt->node_pt(
j) = top_mid_face_node2_pt;
9314 el_pt->node_pt(
j) = top_mid_face_node1_pt;
9323 el_pt->node_pt(
j) = top_mid_face_node0_pt;
9333 for (
int face_index = 0; face_index < 4; face_index++)
9341 face_pt =
new TFace(tet_el_pt->node_pt(1),
9342 tet_el_pt->node_pt(2),
9343 tet_el_pt->node_pt(3));
9348 face_pt =
new TFace(tet_el_pt->node_pt(0),
9349 tet_el_pt->node_pt(2),
9350 tet_el_pt->node_pt(3));
9357 face_pt =
new TFace(tet_el_pt->node_pt(0),
9358 tet_el_pt->node_pt(1),
9359 tet_el_pt->node_pt(3));
9365 face_pt =
new TFace(tet_el_pt->node_pt(0),
9366 tet_el_pt->node_pt(1),
9367 tet_el_pt->node_pt(2));
9372 if (face_pt->is_boundary_face())
9374 std::set<unsigned> bnd;
9375 std::set<unsigned>* bnd_pt = &bnd;
9376 face_pt->get_boundaries_pt(bnd_pt);
9379 if ((*bnd_pt).size() > 1)
9381 std::ostringstream error_stream;
9382 error_stream <<
"TFace should only be on one boundary.\n";
9383 throw OomphLibError(error_stream.str(),
9389 if ((*bnd_pt).size() == 1)
9392 FaceElement* face_el_pt = 0;
9393 if (tet_mesh_is_solid_mesh)
9396 if (
dynamic_cast<SolidTElement<3, 3>*
>(tet_el_pt) == 0)
9398 std::ostringstream error_stream;
9400 <<
"Tet-element cannot be cast to SolidTElement<3,3>.\n"
9401 <<
"BrickFromTetMesh can only be built from\n"
9402 <<
"mesh containing quadratic tets.\n"
9404 throw OomphLibError(error_stream.str(),
9410 face_el_pt =
new DummyFaceElement<SolidTElement<3, 3>>(
9411 tet_el_pt, face_index);
9416 if (
dynamic_cast<TElement<3, 3>*
>(tet_el_pt) == 0)
9418 std::ostringstream error_stream;
9419 error_stream <<
"Tet-element cannot be cast to TElement<3,3>.\n"
9420 <<
"BrickFromTetMesh can only be built from\n"
9421 <<
"mesh containing quadratic tets.\n"
9423 throw OomphLibError(error_stream.str(),
9431 new DummyFaceElement<TElement<3, 3>>(tet_el_pt, face_index);
9437 unsigned b = (*(*bnd_pt).begin());
9439 face_el_pt->set_boundary_number_in_bulk_mesh(
b);
9443 Vector<Node*> brick_face_node_pt(19);
9448 brick_face_node_pt[0] = brick_el1_pt->node_pt(0);
9449 brick_face_node_pt[1] = brick_el3_pt->node_pt(0);
9450 brick_face_node_pt[2] = brick_el2_pt->node_pt(0);
9452 brick_face_node_pt[3] = brick_el1_pt->node_pt(18);
9453 brick_face_node_pt[4] = brick_el2_pt->node_pt(18);
9454 brick_face_node_pt[5] = brick_el1_pt->node_pt(2);
9456 brick_face_node_pt[6] = brick_el1_pt->node_pt(9);
9457 brick_face_node_pt[7] = brick_el3_pt->node_pt(1);
9458 brick_face_node_pt[8] = brick_el3_pt->node_pt(9);
9460 brick_face_node_pt[9] = brick_el2_pt->node_pt(9);
9461 brick_face_node_pt[10] = brick_el2_pt->node_pt(3);
9462 brick_face_node_pt[11] = brick_el1_pt->node_pt(1);
9464 brick_face_node_pt[12] = brick_el1_pt->node_pt(20);
9466 brick_face_node_pt[13] = brick_el2_pt->node_pt(12);
9467 brick_face_node_pt[14] = brick_el1_pt->node_pt(19);
9469 brick_face_node_pt[15] = brick_el1_pt->node_pt(10);
9470 brick_face_node_pt[16] = brick_el2_pt->node_pt(21);
9472 brick_face_node_pt[17] = brick_el3_pt->node_pt(10);
9473 brick_face_node_pt[18] = brick_el1_pt->node_pt(11);
9477 brick_face_node_pt[0] = brick_el0_pt->node_pt(0);
9478 brick_face_node_pt[1] = brick_el3_pt->node_pt(0);
9479 brick_face_node_pt[2] = brick_el2_pt->node_pt(0);
9481 brick_face_node_pt[3] = brick_el0_pt->node_pt(18);
9482 brick_face_node_pt[4] = brick_el2_pt->node_pt(18);
9483 brick_face_node_pt[5] = brick_el0_pt->node_pt(6);
9485 brick_face_node_pt[6] = brick_el0_pt->node_pt(9);
9486 brick_face_node_pt[7] = brick_el3_pt->node_pt(3);
9487 brick_face_node_pt[8] = brick_el3_pt->node_pt(9);
9489 brick_face_node_pt[9] = brick_el2_pt->node_pt(9);
9490 brick_face_node_pt[10] = brick_el2_pt->node_pt(1);
9491 brick_face_node_pt[11] = brick_el0_pt->node_pt(3);
9493 brick_face_node_pt[12] = brick_el0_pt->node_pt(24);
9495 brick_face_node_pt[13] = brick_el2_pt->node_pt(10);
9496 brick_face_node_pt[14] = brick_el0_pt->node_pt(21);
9498 brick_face_node_pt[15] = brick_el0_pt->node_pt(12);
9499 brick_face_node_pt[16] = brick_el3_pt->node_pt(21);
9501 brick_face_node_pt[17] = brick_el3_pt->node_pt(12);
9502 brick_face_node_pt[18] = brick_el0_pt->node_pt(15);
9506 brick_face_node_pt[0] = brick_el0_pt->node_pt(0);
9507 brick_face_node_pt[1] = brick_el1_pt->node_pt(0);
9508 brick_face_node_pt[2] = brick_el2_pt->node_pt(0);
9510 brick_face_node_pt[3] = brick_el0_pt->node_pt(2);
9511 brick_face_node_pt[4] = brick_el1_pt->node_pt(2);
9512 brick_face_node_pt[5] = brick_el0_pt->node_pt(6);
9514 brick_face_node_pt[6] = brick_el0_pt->node_pt(1);
9515 brick_face_node_pt[7] = brick_el1_pt->node_pt(3);
9516 brick_face_node_pt[8] = brick_el1_pt->node_pt(1);
9518 brick_face_node_pt[9] = brick_el2_pt->node_pt(3);
9519 brick_face_node_pt[10] = brick_el2_pt->node_pt(1);
9520 brick_face_node_pt[11] = brick_el0_pt->node_pt(3);
9522 brick_face_node_pt[12] = brick_el0_pt->node_pt(8);
9524 brick_face_node_pt[13] = brick_el2_pt->node_pt(4);
9525 brick_face_node_pt[14] = brick_el0_pt->node_pt(5);
9527 brick_face_node_pt[15] = brick_el0_pt->node_pt(4);
9528 brick_face_node_pt[16] = brick_el1_pt->node_pt(5);
9530 brick_face_node_pt[17] = brick_el1_pt->node_pt(4);
9531 brick_face_node_pt[18] = brick_el0_pt->node_pt(7);
9535 brick_face_node_pt[0] = brick_el1_pt->node_pt(0);
9536 brick_face_node_pt[1] = brick_el3_pt->node_pt(0);
9537 brick_face_node_pt[2] = brick_el0_pt->node_pt(0);
9539 brick_face_node_pt[3] = brick_el1_pt->node_pt(18);
9540 brick_face_node_pt[4] = brick_el3_pt->node_pt(6);
9541 brick_face_node_pt[5] = brick_el1_pt->node_pt(6);
9543 brick_face_node_pt[6] = brick_el1_pt->node_pt(9);
9544 brick_face_node_pt[7] = brick_el3_pt->node_pt(1);
9545 brick_face_node_pt[8] = brick_el3_pt->node_pt(3);
9547 brick_face_node_pt[9] = brick_el0_pt->node_pt(9);
9548 brick_face_node_pt[10] = brick_el0_pt->node_pt(1);
9549 brick_face_node_pt[11] = brick_el1_pt->node_pt(3);
9551 brick_face_node_pt[12] = brick_el1_pt->node_pt(24);
9553 brick_face_node_pt[13] = brick_el0_pt->node_pt(10);
9554 brick_face_node_pt[14] = brick_el1_pt->node_pt(21);
9556 brick_face_node_pt[15] = brick_el1_pt->node_pt(12);
9557 brick_face_node_pt[16] = brick_el3_pt->node_pt(7);
9559 brick_face_node_pt[17] = brick_el3_pt->node_pt(4);
9560 brick_face_node_pt[18] = brick_el1_pt->node_pt(15);
9566 Vector<unsigned> translate(19);
9569 for (
unsigned i = 0;
i < 19;
i++)
9575 for (
unsigned j = 0;
j < 19;
j++)
9578 Node* brick_node_pt = brick_face_node_pt[translate[
j]];
9581 Vector<double>
s = s_face[
j];
9582 Vector<double>
zeta(2);
9583 Vector<double>
x(3);
9584 face_el_pt->interpolated_zeta(
s,
zeta);
9585 face_el_pt->interpolated_x(
s,
x);
9589 double dist =
sqrt(
pow(brick_node_pt->x(0) -
x[0], 2) +
9590 pow(brick_node_pt->x(1) -
x[1], 2) +
9591 pow(brick_node_pt->x(2) -
x[2], 2));
9594 std::ofstream brick0;
9595 std::ofstream brick1;
9596 std::ofstream brick2;
9597 std::ofstream brick3;
9598 brick0.open(
"full_brick0.dat");
9599 brick1.open(
"full_brick1.dat");
9600 brick2.open(
"full_brick2.dat");
9601 brick3.open(
"full_brick3.dat");
9602 for (
unsigned j = 0;
j < 27;
j++)
9604 brick0 << brick_el0_pt->node_pt(
j)->x(0) <<
" "
9605 << brick_el0_pt->node_pt(
j)->x(1) <<
" "
9606 << brick_el0_pt->node_pt(
j)->x(2) <<
"\n";
9608 brick1 << brick_el1_pt->node_pt(
j)->x(0) <<
" "
9609 << brick_el1_pt->node_pt(
j)->x(1) <<
" "
9610 << brick_el1_pt->node_pt(
j)->x(2) <<
"\n";
9612 brick2 << brick_el2_pt->node_pt(
j)->x(0) <<
" "
9613 << brick_el2_pt->node_pt(
j)->x(1) <<
" "
9614 << brick_el2_pt->node_pt(
j)->x(2) <<
"\n";
9616 brick3 << brick_el3_pt->node_pt(
j)->x(0) <<
" "
9617 << brick_el3_pt->node_pt(
j)->x(1) <<
" "
9618 << brick_el3_pt->node_pt(
j)->x(2) <<
"\n";
9625 std::ofstream full_face;
9626 full_face.open(
"full_face.dat");
9627 for (
unsigned j = 0;
j < 6;
j++)
9629 full_face << face_el_pt->node_pt(
j)->x(0) <<
" "
9630 << face_el_pt->node_pt(
j)->x(1) <<
" "
9631 << face_el_pt->node_pt(
j)->x(2) <<
"\n";
9636 int normal_sign = face_el_pt->normal_sign();
9638 std::ostringstream error_stream;
9640 <<
"During assignment of boundary cordinates, the distance\n"
9641 <<
"between brick node and reference point in \n"
9642 <<
"triangular FaceElement is " << dist <<
" which \n"
9643 <<
"is bigger than the tolerance defined in \n"
9644 <<
"BrickFromTetMeshHelper::Face_position_tolerance="
9646 <<
"If this is tolerable, increase the tolerance \n"
9647 <<
"(it's defined in a namespace and therefore publically\n"
9648 <<
"accessible). If not, the Face may be inverted in which \n"
9649 <<
"case you should re-implement the translation scheme,\n"
9650 <<
"following the pattern used in the "
9651 "ThinLayerBrickOnTetMesh."
9652 <<
"\nThe required code fragements are already here but \n"
9653 <<
"the translation is the unit map.\n"
9654 <<
"To aid the diagnostics, the files full_brick[0-3].dat\n"
9655 <<
"contain the coordinates of the 27 nodes in the four\n"
9656 <<
"bricks associated with the current tet and "
9658 <<
"contains the coordinates of the 6 nodes in the "
9660 <<
"\nfrom which the boundary coordinates are extracted.\n"
9661 <<
"FYI: The normal_sign of the face is: " << normal_sign
9663 throw OomphLibError(error_stream.str(),
9671 brick_node_pt->set_coordinates_on_boundary(
b,
zeta);
9739 for (
unsigned e = 0;
e < 4;
e++)
9741 for (
unsigned j = 0;
j < 8;
j++)
9743 delete dummy_q_el_pt[
e]->node_pt(
j);
9745 delete dummy_q_el_pt[
e];
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Scalar * b
Definition: benchVecAdd.cpp:17
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
Vector< Vector< FiniteElement * > > Boundary_element_pt
Definition: mesh.h:91
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:87
std::vector< bool > Boundary_coordinate_exists
Definition: mesh.h:190
Vector< Vector< int > > Face_index_at_boundary
Definition: mesh.h:95
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
RealScalar s
Definition: level1_cplx_impl.h:130
int nb
Definition: level2_impl.h:286
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
double Face_position_tolerance
Tolerance for mismatch during setup of boundary coordinates.
Definition: brick_mesh.cc:37
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 1.0.
Definition: Qelement_face_coordinate_translation_schemes.cc:44
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = -1.0.
Definition: Qelement_face_coordinate_translation_schemes.cc:38
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the east face (s0 = 1.0)
Definition: Qelement_face_coordinate_translation_schemes.cc:93
list x
Definition: plotDoE.py:28
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2