thin_layer_brick_on_tet_mesh.template.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2026 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_THIN_LAYER_BRICK_ON_TET_MESH_TEMPLATE_HEADER
27#define OOMPH_THIN_LAYER_BRICK_ON_TET_MESH_TEMPLATE_HEADER
28
29#ifndef OOMPH_THIN_LAYER_BRICK_ON_TET_MESH_HEADER
30#error __FILE__ should only be included from thin_layer_brick_on_tet_mesh.h.
31#endif // OOMPH_THIN_LAYER_BRICK_ON_TET_MESH_HEADER
32
33#include "../solid/solid_elements.h"
34
35namespace oomph
36{
37 //=====================================================================
38 /// Constructor: Specify (quadratic) tet mesh, boundary IDs of
39 /// boundary on which the current mesh is to be erected (in an FSI context
40 /// this boundary tends to be the FSI boundary of the fluid mesh. Also
41 /// specify the uniform thickness of layer, and the number of element layers.
42 /// The vectors stored in in_out_boundary_ids contain the boundary
43 /// IDs of the other boundaries in the tet mesh. In an FSI context
44 /// these typically identify the in/outflow boundaries in the fluid
45 /// mesh. The boundary enumeration of the current mesh follows the
46 /// one of the underlying fluid mesh: The enumeration of the FSI boundary
47 /// matches (to enable the setup of the FSI matching); the "in/outflow"
48 /// faces in this mesh inherit the same enumeration as the in/outflow
49 /// faces in the underlying fluid mesh. Finally, the "outer" boundary
50 /// gets its own boundary ID.
51 /// Timestepper defaults to steady pseudo-timestepper.
52 //=====================================================================
53 template<class ELEMENT>
54 ThinLayerBrickOnTetMesh<ELEMENT>::ThinLayerBrickOnTetMesh(
55 Mesh* tet_mesh_pt,
56 const Vector<unsigned>& boundary_ids,
57 ThicknessFctPt thickness_fct_pt,
58 const unsigned& nlayer,
59 const Vector<Vector<unsigned>>& in_out_boundary_id,
60 TimeStepper* time_stepper_pt)
61 : Thickness_fct_pt(thickness_fct_pt)
62 {
63 // Mesh can only be built with 3D Qelements.
64 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
65
66 // Figure out if the tet mesh is a solid mesh
67 bool tet_mesh_is_solid_mesh = false;
68 if (dynamic_cast<SolidFiniteElement*>(tet_mesh_pt->element_pt(0)) != 0)
69 {
70 tet_mesh_is_solid_mesh = true;
71 }
72
73 // Setup lookup scheme for local coordinates on triangular faces.
74 // The local coordinates identify the points on the triangular
75 // FaceElements on which we place the bottom layer of the
76 // brick nodes.
77 Vector<Vector<double>> s_face(19);
78 for (unsigned i = 0; i < 19; i++)
79 {
80 s_face[i].resize(2);
81
82 switch (i)
83 {
84 // Vertex nodes
85
86 case 0:
87 s_face[i][0] = 1.0;
88 s_face[i][1] = 0.0;
89 break;
90
91 case 1:
92 s_face[i][0] = 0.0;
93 s_face[i][1] = 1.0;
94 break;
95
96 case 2:
97 s_face[i][0] = 0.0;
98 s_face[i][1] = 0.0;
99 break;
100
101 // Midside nodes
102
103 case 3:
104 s_face[i][0] = 0.5;
105 s_face[i][1] = 0.5;
106 break;
107
108 case 4:
109 s_face[i][0] = 0.0;
110 s_face[i][1] = 0.5;
111 break;
112
113 case 5:
114 s_face[i][0] = 0.5;
115 s_face[i][1] = 0.0;
116 break;
117
118 // Quarter side nodes
119
120 case 6:
121 s_face[i][0] = 0.75;
122 s_face[i][1] = 0.25;
123 break;
124
125 case 7:
126 s_face[i][0] = 0.25;
127 s_face[i][1] = 0.75;
128 break;
129
130 case 8:
131 s_face[i][0] = 0.0;
132 s_face[i][1] = 0.75;
133 break;
134
135 case 9:
136 s_face[i][0] = 0.0;
137 s_face[i][1] = 0.25;
138 break;
139
140 case 10:
141 s_face[i][0] = 0.25;
142 s_face[i][1] = 0.0;
143 break;
144
145 case 11:
146 s_face[i][0] = 0.75;
147 s_face[i][1] = 0.0;
148 break;
149
150 // Central node
151
152 case 12:
153 s_face[i][0] = 1.0 / 3.0;
154 s_face[i][1] = 1.0 / 3.0;
155 break;
156
157 // Vertical internal midside nodes connecting 2 and 3
158
159 case 13:
160 s_face[i][0] = 5.0 / 24.0;
161 s_face[i][1] = 5.0 / 24.0;
162 break;
163
164 case 14:
165 s_face[i][0] = 5.0 / 12.0;
166 s_face[i][1] = 5.0 / 12.0;
167 break;
168
169 // Internal midside nodes connecting nodes 0 and 4
170
171 case 15:
172 s_face[i][1] = 5.0 / 24.0;
173 s_face[i][0] = 7.0 / 12.0; // 1.0-2.0*5.0/24.0;
174 break;
175
176 case 16:
177 s_face[i][1] = 5.0 / 12.0;
178 s_face[i][0] = 1.0 / 6.0; // 1.0-2.0*5.0/12.0;
179 break;
180
181 // Internal midside nodes connecting nodes 1 and 5
182
183 case 17:
184 s_face[i][0] = 5.0 / 24.0;
185 s_face[i][1] = 7.0 / 12.0; // 1.0-2.0*5.0/24.0;
186 break;
187
188 case 18:
189 s_face[i][0] = 5.0 / 12.0;
190 s_face[i][1] = 1.0 / 6.0; // 1.0-2.0*5.0/12.0;
191 break;
192 }
193 }
194
195 // Translation scheme for inverted FaceElements
196 MapMatrixMixed<int, unsigned, unsigned> translate;
197
198 // Initialise with identify mapping
199 for (unsigned i = 0; i < 19; i++)
200 {
201 translate(-1, i) = i;
202 translate(1, i) = i;
203 }
204 translate(-1, 6) = 11;
205 translate(-1, 11) = 6;
206 translate(-1, 3) = 5;
207 translate(-1, 5) = 3;
208 translate(-1, 18) = 14;
209 translate(-1, 14) = 18;
210 translate(-1, 7) = 10;
211 translate(-1, 10) = 7;
212 translate(-1, 13) = 17;
213 translate(-1, 17) = 13;
214 translate(-1, 1) = 2;
215 translate(-1, 2) = 1;
216 translate(-1, 9) = 8;
217 translate(-1, 8) = 9;
218
219 // Lookup scheme relating "fluid" nodes to newly created "solid" nodes
220 // (terminology for fsi problem)
221 std::map<Node*, Node*> solid_node_pt;
222
223 // Look up scheme for quarter edge nodes
224 std::map<Edge, Node*> quarter_edge_node;
225
226 // Map to store normal vectors for all surface nodes, labeled
227 // by node on FSI surface
228 std::map<Node*, Vector<Vector<double>>> normals;
229
230 // Map of nodes connected to node on the tet surface, labeled by
231 // node on FSI surface
232 std::map<Node*, Vector<Node*>> connected_node_pt;
233
234 // Number of elements in brick mesh
235 Element_pt.reserve(3 * boundary_ids.size() * nlayer);
236
237 // Get total number of distinct boundary IDs that we touch
238 // in the fluid mesh
239 std::set<unsigned> all_bnd;
240
241 // Loop over all boundaries in tet mesh that make up the FSI interface
242 unsigned nb = boundary_ids.size();
243 for (unsigned ib = 0; ib < nb; ib++)
244 {
245 // Boundary number in "fluid" tet mesh
246 unsigned b = boundary_ids[ib];
247
248 // Loop over boundary nodes in the fluid mesh on that
249 // boundary
250 unsigned nnod = tet_mesh_pt->nboundary_node(b);
251 for (unsigned j = 0; j < nnod; j++)
252 {
253 Node* nod_pt = tet_mesh_pt->boundary_node_pt(b, j);
254
255 // Get pointer to set of boundaries this node is located on
256 std::set<unsigned>* bnd_pt;
257 nod_pt->get_boundaries_pt(bnd_pt);
258
259 // Add
260 for (std::set<unsigned>::iterator it = (*bnd_pt).begin();
261 it != (*bnd_pt).end();
262 it++)
263 {
264 all_bnd.insert(*it);
265 }
266 }
267 }
268
269 // Highest boundary ID
270 unsigned highest_fluid_bound_id =
271 *std::max_element(all_bnd.begin(), all_bnd.end());
272
273 // Figure out which boundaries are actually on fsi boundary
274 std::vector<bool> is_on_fsi_boundary(highest_fluid_bound_id + 1, false);
275 for (unsigned ib = 0; ib < nb; ib++)
276 {
277 is_on_fsi_boundary[boundary_ids[ib]] = true;
278 }
279
280 // Figure out which boundaries are on the identified in/outflow boundaries
281 unsigned n = in_out_boundary_id.size();
282 Vector<std::vector<bool>> is_on_in_out_boundary(n);
283 Vector<std::set<unsigned>> in_out_boundary_id_set(n);
284 for (unsigned j = 0; j < n; j++)
285 {
286 is_on_in_out_boundary[j].resize(highest_fluid_bound_id + 1, false);
287 unsigned nb = in_out_boundary_id[j].size();
288 for (unsigned ib = 0; ib < nb; ib++)
289 {
290 is_on_in_out_boundary[j][in_out_boundary_id[j][ib]] = true;
291 }
292 }
293
294 // Total number of boundaries: All boundaries that we touch
295 // in the fluid mesh (the FSI boundary and the boundaries
296 // on the in/outflow faces -- we flip these up and use
297 // them for all the boundary faces in adjacent stacks
298 // of solid elements) plus one additional boundary for
299 // the outer boundary.
300 unsigned maxb = highest_fluid_bound_id + 2;
301
302 // Set number of boundaries
303 set_nboundary(maxb);
304
305 // Get ready for boundary lookup scheme
306 Boundary_element_pt.resize(maxb);
307 Face_index_at_boundary.resize(maxb);
308
309 // Loop over all boundaries in tet mesh that make up the FSI interface
310 nb = boundary_ids.size();
311 for (unsigned ib = 0; ib < nb; ib++)
312 {
313 // Boundary number in "fluid" tet mesh
314 unsigned b = boundary_ids[ib];
315
316 // We'll setup boundary coordinates for this one
317 set_boundary_coordinate_exists(b);
318
319 // Remember for future reference
320 FSI_boundary_id.push_back(b);
321
322 // Loop over all elements on this boundary
323 unsigned nel = tet_mesh_pt->nboundary_element(b);
324 for (unsigned e = 0; e < nel; e++)
325 {
326 // Get pointer to the bulk fluid element that is adjacent to boundary b
327 FiniteElement* bulk_elem_pt = tet_mesh_pt->boundary_element_pt(b, e);
328
329 // Find the index of the face of element e along boundary b
330 int face_index = tet_mesh_pt->face_index_at_boundary(b, e);
331
332 // Create new face element
333 FaceElement* face_el_pt = 0;
334 if (tet_mesh_is_solid_mesh)
335 {
336#ifdef PARANOID
337 if (dynamic_cast<SolidTElement<3, 3>*>(bulk_elem_pt) == 0)
338 {
339 std::ostringstream error_stream;
340 error_stream
341 << "Tet-element cannot be cast to SolidTElement<3,3>.\n"
342 << "ThinBrickOnTetMesh can only be erected on mesh containing\n"
343 << "quadratic tets." << std::endl;
344 throw OomphLibError(error_stream.str(),
345 OOMPH_CURRENT_FUNCTION,
346 OOMPH_EXCEPTION_LOCATION);
347 }
348#endif
349 face_el_pt =
350 new DummyFaceElement<SolidTElement<3, 3>>(bulk_elem_pt, face_index);
351 }
352 else
353 {
354#ifdef PARANOID
355 if (dynamic_cast<TElement<3, 3>*>(bulk_elem_pt) == 0)
356 {
357 std::ostringstream error_stream;
358 error_stream
359 << "Tet-element cannot be cast to TElement<3,3>.\n"
360 << "ThinBrickOnTetMesh can only be erected on mesh containing\n"
361 << "quadratic tets." << std::endl;
362 throw OomphLibError(error_stream.str(),
363 OOMPH_CURRENT_FUNCTION,
364 OOMPH_EXCEPTION_LOCATION);
365 }
366#endif
367 face_el_pt =
368 new DummyFaceElement<TElement<3, 3>>(bulk_elem_pt, face_index);
369 }
370
371 // Specify boundary id in bulk mesh (needed to extract
372 // boundary coordinate)
373 face_el_pt->set_boundary_number_in_bulk_mesh(b);
374
375 // Create storage for stack of brick elements
376 Vector<Vector<FiniteElement*>> new_el_pt(3);
377
378 // Sign of normal to detect inversion of FaceElement
379 int normal_sign;
380
381 // Loop over the three bricks that are built on the current
382 //---------------------------------------------------------
383 // triangular face
384 //----------------
385 for (unsigned j = 0; j < 3; j++)
386 {
387 // Build stack of bricks
388 new_el_pt[j].resize(nlayer);
389 for (unsigned ilayer = 0; ilayer < nlayer; ilayer++)
390 {
391 new_el_pt[j][ilayer] = new ELEMENT;
392 Element_pt.push_back(new_el_pt[j][ilayer]);
393 }
394
395 Boundary_element_pt[b].push_back(new_el_pt[j][0]);
396 Face_index_at_boundary[b].push_back(-3);
397
398 // Associate zero-th node with vertex of triangular face
399 //------------------------------------------------------
400 unsigned j_local = 0;
401
402 // Get normal sign
403 normal_sign = face_el_pt->normal_sign();
404
405 // Get coordinates etc of point from face: Vertex nodes enumerated
406 // first....
407 Vector<double> s = s_face[translate(normal_sign, j)];
408 Vector<double> zeta(2);
409 Vector<double> x(3);
410 Vector<double> unit_normal(3);
411 face_el_pt->interpolated_zeta(s, zeta);
412 face_el_pt->interpolated_x(s, x);
413 face_el_pt->outer_unit_normal(s, unit_normal);
414
415 // Get node in the "fluid" mesh from face
416 Node* fluid_node_pt = face_el_pt->node_pt(translate(normal_sign, j));
417
418 // Has the corresponding "solid" node already been created?
419 Node* existing_node_pt = solid_node_pt[fluid_node_pt];
420 if (existing_node_pt == 0)
421 {
422 // Create new node
423 Node* new_node_pt = new_el_pt[j][0]->construct_boundary_node(
424 j_local, time_stepper_pt);
425 Node_pt.push_back(new_node_pt);
426
427 //...and remember it
428 solid_node_pt[fluid_node_pt] = new_node_pt;
429
430 // Set coordinates
431 new_node_pt->x(0) = x[0];
432 new_node_pt->x(1) = x[1];
433 new_node_pt->x(2) = x[2];
434
435 // Set boundary stuff -- boundary IDs copied from fluid
436 bool only_on_fsi = true;
437 std::set<unsigned>* bnd_pt;
438 fluid_node_pt->get_boundaries_pt(bnd_pt);
439 for (std::set<unsigned>::iterator it = (*bnd_pt).begin();
440 it != (*bnd_pt).end();
441 it++)
442 {
443 if (!is_on_fsi_boundary[(*it)]) only_on_fsi = false;
444 add_boundary_node((*it), new_node_pt);
445 }
446 new_node_pt->set_coordinates_on_boundary(b, zeta);
447 normals[new_node_pt].push_back(unit_normal);
448
449 // If bottom node is only on FSI boundary, the nodes above
450 // are not boundary nodes, apart from the last one!
451 if (only_on_fsi)
452 {
453 // Create other nodes in bottom layer
454 Node* new_nod_pt =
455 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
456 connected_node_pt[new_node_pt].push_back(new_nod_pt);
457 Node_pt.push_back(new_nod_pt);
458
459 // One layer thick?
460 if (nlayer == 1)
461 {
462 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
463 j_local + 18, time_stepper_pt);
464 connected_node_pt[new_node_pt].push_back(new_nod_pt);
465 Node_pt.push_back(new_nod_pt);
466 }
467 else
468 {
469 Node* new_nod_pt = new_el_pt[j][0]->construct_node(
470 j_local + 18, time_stepper_pt);
471 connected_node_pt[new_node_pt].push_back(new_nod_pt);
472 Node_pt.push_back(new_nod_pt);
473 }
474
475 // Now do other layers
476 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
477 {
478 // Copy bottom node from below
479 new_el_pt[j][ilayer]->node_pt(j_local) =
480 connected_node_pt[new_node_pt][2 * ilayer - 1];
481
482 // Create new nodes
483 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
484 j_local + 9, time_stepper_pt);
485 connected_node_pt[new_node_pt].push_back(new_nod_pt);
486 Node_pt.push_back(new_nod_pt);
487
488 // Last node is boundary node
489 if (ilayer != (nlayer - 1))
490 {
491 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
492 j_local + 18, time_stepper_pt);
493 connected_node_pt[new_node_pt].push_back(new_nod_pt);
494 Node_pt.push_back(new_nod_pt);
495 }
496 else
497 {
498 Node* new_nod_pt =
499 new_el_pt[j][ilayer]->construct_boundary_node(
500 j_local + 18, time_stepper_pt);
501 connected_node_pt[new_node_pt].push_back(new_nod_pt);
502 Node_pt.push_back(new_nod_pt);
503 }
504 }
505 }
506 else
507 {
508 // Create other boundary nodes in bottom layer
509 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
510 j_local + 9, time_stepper_pt);
511 connected_node_pt[new_node_pt].push_back(new_nod_pt);
512 Node_pt.push_back(new_nod_pt);
513
514 new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
515 j_local + 18, time_stepper_pt);
516 connected_node_pt[new_node_pt].push_back(new_nod_pt);
517 Node_pt.push_back(new_nod_pt);
518
519 // Now do other layers
520 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
521 {
522 // Copy bottom node from below
523 new_el_pt[j][ilayer]->node_pt(j_local) =
524 connected_node_pt[new_node_pt][2 * ilayer - 1];
525
526 // Create new boundary nodes
527 Node* new_nod_pt =
528 new_el_pt[j][ilayer]->construct_boundary_node(
529 j_local + 9, time_stepper_pt);
530 connected_node_pt[new_node_pt].push_back(new_nod_pt);
531 Node_pt.push_back(new_nod_pt);
532 new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
533 j_local + 18, time_stepper_pt);
534 connected_node_pt[new_node_pt].push_back(new_nod_pt);
535 Node_pt.push_back(new_nod_pt);
536 }
537 }
538 }
539 else
540 {
541 // Add (repeated) bottom node to its other boundary and add
542 // coordinates
543 existing_node_pt->set_coordinates_on_boundary(b, zeta);
544 normals[existing_node_pt].push_back(unit_normal);
545
546 // Get pointer to nodes in bottom layer
547 new_el_pt[j][0]->node_pt(j_local) = existing_node_pt;
548 new_el_pt[j][0]->node_pt(j_local + 9) =
549 connected_node_pt[existing_node_pt][0];
550 new_el_pt[j][0]->node_pt(j_local + 18) =
551 connected_node_pt[existing_node_pt][1];
552
553 // Now do other layers
554 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
555 {
556 new_el_pt[j][ilayer]->node_pt(j_local) =
557 connected_node_pt[existing_node_pt][2 * ilayer - 1];
558 new_el_pt[j][ilayer]->node_pt(j_local + 9) =
559 connected_node_pt[existing_node_pt][2 * ilayer];
560 new_el_pt[j][ilayer]->node_pt(j_local + 18) =
561 connected_node_pt[existing_node_pt][2 * ilayer + 1];
562 }
563 }
564
565 // Second node with midside node in triangular face
566 //-------------------------------------------------
567 j_local = 2;
568
569 // Get coordinates etc of point from face: Midside nodes enumerated
570 // after vertex nodes
571 s = s_face[translate(normal_sign, j + 3)];
572 face_el_pt->interpolated_zeta(s, zeta);
573 face_el_pt->interpolated_x(s, x);
574 face_el_pt->outer_unit_normal(s, unit_normal);
575
576 // Get node in the "fluid" mesh from face
577 fluid_node_pt = face_el_pt->node_pt(translate(normal_sign, j + 3));
578
579 // Has the corresponding "solid" node already been created?
580 existing_node_pt = solid_node_pt[fluid_node_pt];
581 if (existing_node_pt == 0)
582 {
583 // Create new node
584 Node* new_node_pt = new_el_pt[j][0]->construct_boundary_node(
585 j_local, time_stepper_pt);
586 Node_pt.push_back(new_node_pt);
587
588 // ...and remember it
589 solid_node_pt[fluid_node_pt] = new_node_pt;
590
591 // Set coordinates
592 new_node_pt->x(0) = x[0];
593 new_node_pt->x(1) = x[1];
594 new_node_pt->x(2) = x[2];
595
596 // Set boundary stuff -- boundary IDs copied from fluid
597 bool only_on_fsi = true;
598 std::set<unsigned>* bnd_pt;
599 fluid_node_pt->get_boundaries_pt(bnd_pt);
600 for (std::set<unsigned>::iterator it = (*bnd_pt).begin();
601 it != (*bnd_pt).end();
602 it++)
603 {
604 if (!is_on_fsi_boundary[(*it)]) only_on_fsi = false;
605 add_boundary_node((*it), new_node_pt);
606 }
607 new_node_pt->set_coordinates_on_boundary(b, zeta);
608 normals[new_node_pt].push_back(unit_normal);
609
610 // If bottom node is only on FSI boundary, the nodes above
611 // are not boundary nodes, apart from the last one!
612 if (only_on_fsi)
613 {
614 // Create other nodes in bottom layer
615 Node* new_nod_pt =
616 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
617 connected_node_pt[new_node_pt].push_back(new_nod_pt);
618 Node_pt.push_back(new_nod_pt);
619
620 // One layer thick?
621 if (nlayer == 1)
622 {
623 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
624 j_local + 18, time_stepper_pt);
625 connected_node_pt[new_node_pt].push_back(new_nod_pt);
626 Node_pt.push_back(new_nod_pt);
627 }
628 else
629 {
630 Node* new_nod_pt = new_el_pt[j][0]->construct_node(
631 j_local + 18, time_stepper_pt);
632 connected_node_pt[new_node_pt].push_back(new_nod_pt);
633 Node_pt.push_back(new_nod_pt);
634 }
635
636 // Now do other layers
637 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
638 {
639 // Copy bottom node from below
640 new_el_pt[j][ilayer]->node_pt(j_local) =
641 connected_node_pt[new_node_pt][2 * ilayer - 1];
642
643 // Create new nodes
644 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
645 j_local + 9, time_stepper_pt);
646 connected_node_pt[new_node_pt].push_back(new_nod_pt);
647 Node_pt.push_back(new_nod_pt);
648
649 // Last node is boundary node
650 if (ilayer != (nlayer - 1))
651 {
652 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
653 j_local + 18, time_stepper_pt);
654 connected_node_pt[new_node_pt].push_back(new_nod_pt);
655 Node_pt.push_back(new_nod_pt);
656 }
657 else
658 {
659 Node* new_nod_pt =
660 new_el_pt[j][ilayer]->construct_boundary_node(
661 j_local + 18, time_stepper_pt);
662 connected_node_pt[new_node_pt].push_back(new_nod_pt);
663 Node_pt.push_back(new_nod_pt);
664 }
665 }
666 }
667 else
668 {
669 // Create other boundary nodes in bottom layer
670 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
671 j_local + 9, time_stepper_pt);
672 connected_node_pt[new_node_pt].push_back(new_nod_pt);
673 Node_pt.push_back(new_nod_pt);
674
675 new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
676 j_local + 18, time_stepper_pt);
677 connected_node_pt[new_node_pt].push_back(new_nod_pt);
678 Node_pt.push_back(new_nod_pt);
679
680 // Now do other layers
681 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
682 {
683 // Copy bottom node from below
684 new_el_pt[j][ilayer]->node_pt(j_local) =
685 connected_node_pt[new_node_pt][2 * ilayer - 1];
686
687 // Create new boundary nodes
688 Node* new_nod_pt =
689 new_el_pt[j][ilayer]->construct_boundary_node(
690 j_local + 9, time_stepper_pt);
691 connected_node_pt[new_node_pt].push_back(new_nod_pt);
692 Node_pt.push_back(new_nod_pt);
693 new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
694 j_local + 18, time_stepper_pt);
695 connected_node_pt[new_node_pt].push_back(new_nod_pt);
696 Node_pt.push_back(new_nod_pt);
697 }
698 }
699 }
700 else
701 {
702 // Add (repeated) bottom node to its other boundary and add
703 // coordinates
704 existing_node_pt->set_coordinates_on_boundary(b, zeta);
705 normals[existing_node_pt].push_back(unit_normal);
706
707 // Get pointer to nodes in bottom layer
708 new_el_pt[j][0]->node_pt(j_local) = existing_node_pt;
709 new_el_pt[j][0]->node_pt(j_local + 9) =
710 connected_node_pt[existing_node_pt][0];
711 new_el_pt[j][0]->node_pt(j_local + 18) =
712 connected_node_pt[existing_node_pt][1];
713
714 // Now do other layers
715 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
716 {
717 new_el_pt[j][ilayer]->node_pt(j_local) =
718 connected_node_pt[existing_node_pt][2 * ilayer - 1];
719 new_el_pt[j][ilayer]->node_pt(j_local + 9) =
720 connected_node_pt[existing_node_pt][2 * ilayer];
721 new_el_pt[j][ilayer]->node_pt(j_local + 18) =
722 connected_node_pt[existing_node_pt][2 * ilayer + 1];
723 }
724 }
725
726 // First node is quarter-edge node on triangular face
727 //---------------------------------------------------
728 j_local = 1;
729
730 // Get coordinates of point from face: Quarter edge nodes enumerated
731 // after midside nodes
732 s = s_face[translate(normal_sign, 6 + 2 * j)];
733 face_el_pt->interpolated_zeta(s, zeta);
734 face_el_pt->interpolated_x(s, x);
735 face_el_pt->outer_unit_normal(s, unit_normal);
736
737 // Create Edge
738 Edge edge(face_el_pt->node_pt(translate(normal_sign, j)),
739 face_el_pt->node_pt(translate(normal_sign, j + 3)));
740
741 // Does node already exist?
742 existing_node_pt = quarter_edge_node[edge];
743 if (existing_node_pt == 0)
744 {
745 // Create new node
746 Node* new_node_pt = new_el_pt[j][0]->construct_boundary_node(
747 j_local, time_stepper_pt);
748 Node_pt.push_back(new_node_pt);
749
750 //...and remember it
751 quarter_edge_node[edge] = new_node_pt;
752
753 // Set coordinates
754 new_node_pt->x(0) = x[0];
755 new_node_pt->x(1) = x[1];
756 new_node_pt->x(2) = x[2];
757
758 // Set boundary stuff -- boundary IDs copied from fluid
759 std::set<unsigned>* bnd1_pt;
760 edge.node1_pt()->get_boundaries_pt(bnd1_pt);
761 std::set<unsigned>* bnd2_pt;
762 edge.node2_pt()->get_boundaries_pt(bnd2_pt);
763 std::set<unsigned> bnd;
764 set_intersection((*bnd1_pt).begin(),
765 (*bnd1_pt).end(),
766 (*bnd2_pt).begin(),
767 (*bnd2_pt).end(),
768 inserter(bnd, bnd.begin()));
769 bool only_on_fsi = true;
770 for (std::set<unsigned>::iterator it = bnd.begin(); it != bnd.end();
771 it++)
772 {
773 if (!is_on_fsi_boundary[(*it)]) only_on_fsi = false;
774 add_boundary_node((*it), new_node_pt);
775 }
776 new_node_pt->set_coordinates_on_boundary(b, zeta);
777 normals[new_node_pt].push_back(unit_normal);
778
779 // If bottom node is only on FSI boundary, the nodes above
780 // are not boundary nodes, apart from the last one!
781 if (only_on_fsi)
782 {
783 // Create other nodes in bottom layer
784 Node* new_nod_pt =
785 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
786 connected_node_pt[new_node_pt].push_back(new_nod_pt);
787 Node_pt.push_back(new_nod_pt);
788
789 // One layer thick?
790 if (nlayer == 1)
791 {
792 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
793 j_local + 18, time_stepper_pt);
794 connected_node_pt[new_node_pt].push_back(new_nod_pt);
795 Node_pt.push_back(new_nod_pt);
796 }
797 else
798 {
799 Node* new_nod_pt = new_el_pt[j][0]->construct_node(
800 j_local + 18, time_stepper_pt);
801 connected_node_pt[new_node_pt].push_back(new_nod_pt);
802 Node_pt.push_back(new_nod_pt);
803 }
804
805 // Now do other layers
806 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
807 {
808 // Copy bottom node from below
809 new_el_pt[j][ilayer]->node_pt(j_local) =
810 connected_node_pt[new_node_pt][2 * ilayer - 1];
811
812 // Create new nodes
813 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
814 j_local + 9, time_stepper_pt);
815 connected_node_pt[new_node_pt].push_back(new_nod_pt);
816 Node_pt.push_back(new_nod_pt);
817
818 // Last node is boundary node
819 if (ilayer != (nlayer - 1))
820 {
821 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
822 j_local + 18, time_stepper_pt);
823 connected_node_pt[new_node_pt].push_back(new_nod_pt);
824 Node_pt.push_back(new_nod_pt);
825 }
826 else
827 {
828 Node* new_nod_pt =
829 new_el_pt[j][ilayer]->construct_boundary_node(
830 j_local + 18, time_stepper_pt);
831 connected_node_pt[new_node_pt].push_back(new_nod_pt);
832 Node_pt.push_back(new_nod_pt);
833 }
834 }
835 }
836 else
837 {
838 // Create other boundary nodes in bottom layer
839 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
840 j_local + 9, time_stepper_pt);
841 connected_node_pt[new_node_pt].push_back(new_nod_pt);
842 Node_pt.push_back(new_nod_pt);
843
844 new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
845 j_local + 18, time_stepper_pt);
846 connected_node_pt[new_node_pt].push_back(new_nod_pt);
847 Node_pt.push_back(new_nod_pt);
848
849 // Now do other layers
850 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
851 {
852 // Copy bottom node from below
853 new_el_pt[j][ilayer]->node_pt(j_local) =
854 connected_node_pt[new_node_pt][2 * ilayer - 1];
855
856 // Create new boundary nodes
857 Node* new_nod_pt =
858 new_el_pt[j][ilayer]->construct_boundary_node(
859 j_local + 9, time_stepper_pt);
860 connected_node_pt[new_node_pt].push_back(new_nod_pt);
861 Node_pt.push_back(new_nod_pt);
862 new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
863 j_local + 18, time_stepper_pt);
864 connected_node_pt[new_node_pt].push_back(new_nod_pt);
865 Node_pt.push_back(new_nod_pt);
866 }
867 }
868 }
869 else
870 {
871 // Add (repeated) bottom node to its other boundary and add
872 // coordinates
873 existing_node_pt->set_coordinates_on_boundary(b, zeta);
874 normals[existing_node_pt].push_back(unit_normal);
875
876 // Get pointer to nodes in bottom layer
877 new_el_pt[j][0]->node_pt(j_local) = existing_node_pt;
878 new_el_pt[j][0]->node_pt(j_local + 9) =
879 connected_node_pt[existing_node_pt][0];
880 new_el_pt[j][0]->node_pt(j_local + 18) =
881 connected_node_pt[existing_node_pt][1];
882
883 // Now do other layers
884 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
885 {
886 new_el_pt[j][ilayer]->node_pt(j_local) =
887 connected_node_pt[existing_node_pt][2 * ilayer - 1];
888 new_el_pt[j][ilayer]->node_pt(j_local + 9) =
889 connected_node_pt[existing_node_pt][2 * ilayer];
890 new_el_pt[j][ilayer]->node_pt(j_local + 18) =
891 connected_node_pt[existing_node_pt][2 * ilayer + 1];
892 }
893 }
894
895 // Third node is three-quarter-edge node on triangular face
896 //---------------------------------------------------------
897 j_local = 3;
898
899 // Create Edge
900 unsigned other_node = 0;
901 unsigned jj = 0;
902 switch (j)
903 {
904 case 0:
905 other_node = 5;
906 jj = 11;
907 break;
908 case 1:
909 other_node = 3;
910 jj = 7;
911 break;
912 case 2:
913 other_node = 4;
914 jj = 9;
915 break;
916 }
917 Edge edge2(face_el_pt->node_pt(translate(normal_sign, j)),
918 face_el_pt->node_pt(translate(normal_sign, other_node)));
919
920 // Get coordinates of point from face:
921 s = s_face[translate(normal_sign, jj)];
922 face_el_pt->interpolated_zeta(s, zeta);
923 face_el_pt->interpolated_x(s, x);
924 face_el_pt->outer_unit_normal(s, unit_normal);
925
926 // Does node already exist?
927 existing_node_pt = quarter_edge_node[edge2];
928 if (existing_node_pt == 0)
929 {
930 // Create new node
931 Node* new_node_pt = new_el_pt[j][0]->construct_boundary_node(
932 j_local, time_stepper_pt);
933 Node_pt.push_back(new_node_pt);
934
935 //..and remember it
936 quarter_edge_node[edge2] = new_node_pt;
937
938 // Set coordinates
939 new_node_pt->x(0) = x[0];
940 new_node_pt->x(1) = x[1];
941 new_node_pt->x(2) = x[2];
942
943 // Set boundary stuff -- boundary IDs copied from fluid
944 std::set<unsigned>* bnd1_pt;
945 edge2.node1_pt()->get_boundaries_pt(bnd1_pt);
946 std::set<unsigned>* bnd2_pt;
947 edge2.node2_pt()->get_boundaries_pt(bnd2_pt);
948 std::set<unsigned> bnd;
949 set_intersection((*bnd1_pt).begin(),
950 (*bnd1_pt).end(),
951 (*bnd2_pt).begin(),
952 (*bnd2_pt).end(),
953 inserter(bnd, bnd.begin()));
954 bool only_on_fsi = true;
955 for (std::set<unsigned>::iterator it = bnd.begin(); it != bnd.end();
956 it++)
957 {
958 if (!is_on_fsi_boundary[(*it)]) only_on_fsi = false;
959 add_boundary_node((*it), new_node_pt);
960 }
961 new_node_pt->set_coordinates_on_boundary(b, zeta);
962 normals[new_node_pt].push_back(unit_normal);
963
964 // If bottom node is only on FSI boundary, the nodes above
965 // are not boundary nodes, apart from the last one!
966 if (only_on_fsi)
967 {
968 // Create other nodes in bottom layer
969 Node* new_nod_pt =
970 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
971 connected_node_pt[new_node_pt].push_back(new_nod_pt);
972 Node_pt.push_back(new_nod_pt);
973
974 // One layer thick?
975 if (nlayer == 1)
976 {
977 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
978 j_local + 18, time_stepper_pt);
979 connected_node_pt[new_node_pt].push_back(new_nod_pt);
980 Node_pt.push_back(new_nod_pt);
981 }
982 else
983 {
984 Node* new_nod_pt = new_el_pt[j][0]->construct_node(
985 j_local + 18, time_stepper_pt);
986 connected_node_pt[new_node_pt].push_back(new_nod_pt);
987 Node_pt.push_back(new_nod_pt);
988 }
989
990 // Now do other layers
991 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
992 {
993 // Copy bottom node from below
994 new_el_pt[j][ilayer]->node_pt(j_local) =
995 connected_node_pt[new_node_pt][2 * ilayer - 1];
996
997 // Create new nodes
998 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
999 j_local + 9, time_stepper_pt);
1000 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1001 Node_pt.push_back(new_nod_pt);
1002 // Last node is boundary node
1003 if (ilayer != (nlayer - 1))
1004 {
1005 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1006 j_local + 18, time_stepper_pt);
1007 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1008 Node_pt.push_back(new_nod_pt);
1009 }
1010 else
1011 {
1012 Node* new_nod_pt =
1013 new_el_pt[j][ilayer]->construct_boundary_node(
1014 j_local + 18, time_stepper_pt);
1015 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1016 Node_pt.push_back(new_nod_pt);
1017 }
1018 }
1019 }
1020 else
1021 {
1022 // Create other boundary nodes in bottom layer
1023 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
1024 j_local + 9, time_stepper_pt);
1025 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1026 Node_pt.push_back(new_nod_pt);
1027
1028 new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
1029 j_local + 18, time_stepper_pt);
1030 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1031 Node_pt.push_back(new_nod_pt);
1032
1033 // Now do other layers
1034 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1035 {
1036 // Copy bottom node from below
1037 new_el_pt[j][ilayer]->node_pt(j_local) =
1038 connected_node_pt[new_node_pt][2 * ilayer - 1];
1039
1040 // Create new boundary nodes
1041 Node* new_nod_pt =
1042 new_el_pt[j][ilayer]->construct_boundary_node(
1043 j_local + 9, time_stepper_pt);
1044 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1045 Node_pt.push_back(new_nod_pt);
1046 new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
1047 j_local + 18, time_stepper_pt);
1048 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1049 Node_pt.push_back(new_nod_pt);
1050 }
1051 }
1052 }
1053 else
1054 {
1055 // Add (repeated) bottom node to its other boundary and add
1056 // coordinates
1057 existing_node_pt->set_coordinates_on_boundary(b, zeta);
1058 normals[existing_node_pt].push_back(unit_normal);
1059
1060 // Get pointer to nodes in bottom layer
1061 new_el_pt[j][0]->node_pt(j_local) = existing_node_pt;
1062 new_el_pt[j][0]->node_pt(j_local + 9) =
1063 connected_node_pt[existing_node_pt][0];
1064 new_el_pt[j][0]->node_pt(j_local + 18) =
1065 connected_node_pt[existing_node_pt][1];
1066
1067 // Now do other layers
1068 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1069 {
1070 new_el_pt[j][ilayer]->node_pt(j_local) =
1071 connected_node_pt[existing_node_pt][2 * ilayer - 1];
1072 new_el_pt[j][ilayer]->node_pt(j_local + 9) =
1073 connected_node_pt[existing_node_pt][2 * ilayer];
1074 new_el_pt[j][ilayer]->node_pt(j_local + 18) =
1075 connected_node_pt[existing_node_pt][2 * ilayer + 1];
1076 }
1077 }
1078
1079 // Fourth node is unique for all elements
1080 //--------------------------------------
1081 j_local = 4;
1082
1083 // Create new node
1084 Node* new_node_pt =
1085 new_el_pt[j][0]->construct_boundary_node(j_local, time_stepper_pt);
1086 Node_pt.push_back(new_node_pt);
1087
1088 jj = 0;
1089 switch (j)
1090 {
1091 case 0:
1092 jj = 15;
1093 break;
1094 case 1:
1095 jj = 17;
1096 break;
1097 case 2:
1098 jj = 13;
1099 break;
1100 }
1101
1102 // Get coordinates etc of point from face:
1103 s = s_face[translate(normal_sign, jj)];
1104 face_el_pt->interpolated_zeta(s, zeta);
1105 face_el_pt->interpolated_x(s, x);
1106 face_el_pt->outer_unit_normal(s, unit_normal);
1107
1108 // Set coordinates
1109 new_node_pt->x(0) = x[0];
1110 new_node_pt->x(1) = x[1];
1111 new_node_pt->x(2) = x[2];
1112
1113 // Set boundary stuff
1114 add_boundary_node(b, new_node_pt);
1115 new_node_pt->set_coordinates_on_boundary(b, zeta);
1116 normals[new_node_pt].push_back(unit_normal);
1117
1118 // Create other nodes in bottom layer
1119 Node* new_nod_pt =
1120 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
1121 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1122 Node_pt.push_back(new_nod_pt);
1123
1124 // One layer thick?
1125 if (nlayer == 1)
1126 {
1127 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
1128 j_local + 18, time_stepper_pt);
1129 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1130 Node_pt.push_back(new_nod_pt);
1131 }
1132 else
1133 {
1134 Node* new_nod_pt =
1135 new_el_pt[j][0]->construct_node(j_local + 18, time_stepper_pt);
1136 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1137 Node_pt.push_back(new_nod_pt);
1138 }
1139
1140 // Now do other layers
1141 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1142 {
1143 // Copy bottom node from below
1144 new_el_pt[j][ilayer]->node_pt(j_local) =
1145 connected_node_pt[new_node_pt][2 * ilayer - 1];
1146
1147 // Create new nodes
1148 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1149 j_local + 9, time_stepper_pt);
1150 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1151 Node_pt.push_back(new_nod_pt);
1152 // Last node is boundary node
1153 if (ilayer != (nlayer - 1))
1154 {
1155 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1156 j_local + 18, time_stepper_pt);
1157 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1158 Node_pt.push_back(new_nod_pt);
1159 }
1160 else
1161 {
1162 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
1163 j_local + 18, time_stepper_pt);
1164 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1165 Node_pt.push_back(new_nod_pt);
1166 }
1167 }
1168
1169 // Fifth node is created by all elements (internal to this
1170 //--------------------------------------------------------
1171 // this patch of elements to can't have been created elsewhere)
1172 //-------------------------------------------------------------
1173
1174 j_local = 5;
1175
1176 // Create new node
1177 new_node_pt =
1178 new_el_pt[j][0]->construct_boundary_node(j_local, time_stepper_pt);
1179 Node_pt.push_back(new_node_pt);
1180
1181 // Get coordinates of point from face:
1182 jj = 0;
1183 switch (j)
1184 {
1185 case 0:
1186 jj = 14;
1187 break;
1188 case 1:
1189 jj = 16;
1190 break;
1191 case 2:
1192 jj = 18;
1193 break;
1194 }
1195
1196 // Get coordinates etc from face
1197 s = s_face[translate(normal_sign, jj)];
1198 face_el_pt->interpolated_zeta(s, zeta);
1199 face_el_pt->interpolated_x(s, x);
1200 face_el_pt->outer_unit_normal(s, unit_normal);
1201
1202 // Set coordinates
1203 new_node_pt->x(0) = x[0];
1204 new_node_pt->x(1) = x[1];
1205 new_node_pt->x(2) = x[2];
1206
1207 // Set boundary stuff
1208 add_boundary_node(b, new_node_pt);
1209 new_node_pt->set_coordinates_on_boundary(b, zeta);
1210 normals[new_node_pt].push_back(unit_normal);
1211
1212 // Create other nodes in bottom layer
1213 new_nod_pt =
1214 new_el_pt[j][0]->construct_node(j_local + 9, time_stepper_pt);
1215 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1216 Node_pt.push_back(new_nod_pt);
1217
1218 // One layer thick?
1219 if (nlayer == 1)
1220 {
1221 Node* new_nod_pt = new_el_pt[j][0]->construct_boundary_node(
1222 j_local + 18, time_stepper_pt);
1223 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1224 Node_pt.push_back(new_nod_pt);
1225 }
1226 else
1227 {
1228 Node* new_nod_pt =
1229 new_el_pt[j][0]->construct_node(j_local + 18, time_stepper_pt);
1230 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1231 Node_pt.push_back(new_nod_pt);
1232 }
1233
1234 // Now do other layers
1235 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1236 {
1237 // Copy bottom node from below
1238 new_el_pt[j][ilayer]->node_pt(j_local) =
1239 connected_node_pt[new_node_pt][2 * ilayer - 1];
1240
1241 // Create other nodes
1242 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1243 j_local + 9, time_stepper_pt);
1244 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1245 Node_pt.push_back(new_nod_pt);
1246
1247 // Last node is boundary node
1248 if (ilayer != (nlayer - 1))
1249 {
1250 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_node(
1251 j_local + 18, time_stepper_pt);
1252 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1253 Node_pt.push_back(new_nod_pt);
1254 }
1255 else
1256 {
1257 Node* new_nod_pt = new_el_pt[j][ilayer]->construct_boundary_node(
1258 j_local + 18, time_stepper_pt);
1259 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1260 Node_pt.push_back(new_nod_pt);
1261 }
1262 }
1263
1264 } // End over the three bricks erected on current triangular face
1265
1266 // Last element builds central node as its node 8
1267 //-----------------------------------------------
1268
1269 unsigned j_local = 8;
1270
1271 // Create new node
1272 Node* new_node_pt =
1273 new_el_pt[2][0]->construct_boundary_node(j_local, time_stepper_pt);
1274 Node_pt.push_back(new_node_pt);
1275
1276 // Get coordinates etc of point from face: Central node is
1277 // node 12 in face enumeration.
1278 Vector<double> s = s_face[12];
1279 Vector<double> zeta(2);
1280 Vector<double> x(3);
1281 Vector<double> unit_normal(3);
1282 face_el_pt->interpolated_zeta(s, zeta);
1283 face_el_pt->interpolated_x(s, x);
1284 face_el_pt->outer_unit_normal(s, unit_normal);
1285
1286 // Set coordinates
1287 new_node_pt->x(0) = x[0];
1288 new_node_pt->x(1) = x[1];
1289 new_node_pt->x(2) = x[2];
1290
1291 // Set boundary stuff
1292 add_boundary_node(b, new_node_pt);
1293 new_node_pt->set_coordinates_on_boundary(b, zeta);
1294 normals[new_node_pt].push_back(unit_normal);
1295
1296 // Create other nodes in bottom layer
1297 Node* new_nod_pt =
1298 new_el_pt[2][0]->construct_node(j_local + 9, time_stepper_pt);
1299 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1300 Node_pt.push_back(new_nod_pt);
1301
1302 // One layer thick?
1303 if (nlayer == 1)
1304 {
1305 Node* new_nod_pt = new_el_pt[2][0]->construct_boundary_node(
1306 j_local + 18, time_stepper_pt);
1307 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1308 Node_pt.push_back(new_nod_pt);
1309 }
1310 else
1311 {
1312 Node* new_nod_pt =
1313 new_el_pt[2][0]->construct_node(j_local + 18, time_stepper_pt);
1314 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1315 Node_pt.push_back(new_nod_pt);
1316 }
1317
1318 // Now do other layers
1319 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1320 {
1321 // Copy bottom node from below
1322 new_el_pt[2][ilayer]->node_pt(j_local) =
1323 connected_node_pt[new_node_pt][2 * ilayer - 1];
1324
1325 // Create other nodes
1326 Node* new_nod_pt =
1327 new_el_pt[2][ilayer]->construct_node(j_local + 9, time_stepper_pt);
1328 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1329 Node_pt.push_back(new_nod_pt);
1330
1331 // Last node is boundary node
1332 if (ilayer != (nlayer - 1))
1333 {
1334 Node* new_nod_pt = new_el_pt[2][ilayer]->construct_node(
1335 j_local + 18, time_stepper_pt);
1336 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1337 Node_pt.push_back(new_nod_pt);
1338 }
1339 else
1340 {
1341 Node* new_nod_pt = new_el_pt[2][ilayer]->construct_boundary_node(
1342 j_local + 18, time_stepper_pt);
1343 connected_node_pt[new_node_pt].push_back(new_nod_pt);
1344 Node_pt.push_back(new_nod_pt);
1345 }
1346 }
1347
1348 // Other elements copy that node across
1349 new_el_pt[1][0]->node_pt(j_local) = new_node_pt;
1350 new_el_pt[0][0]->node_pt(j_local) = new_node_pt;
1351
1352 new_el_pt[1][0]->node_pt(j_local + 9) =
1353 connected_node_pt[new_node_pt][0];
1354 new_el_pt[0][0]->node_pt(j_local + 9) =
1355 connected_node_pt[new_node_pt][0];
1356
1357 new_el_pt[1][0]->node_pt(j_local + 18) =
1358 connected_node_pt[new_node_pt][1];
1359 new_el_pt[0][0]->node_pt(j_local + 18) =
1360 connected_node_pt[new_node_pt][1];
1361
1362 // Now do layers
1363 for (unsigned ilayer = 1; ilayer < nlayer; ilayer++)
1364 {
1365 new_el_pt[1][ilayer]->node_pt(j_local) =
1366 connected_node_pt[new_node_pt][2 * ilayer - 1];
1367 new_el_pt[0][ilayer]->node_pt(j_local) =
1368 connected_node_pt[new_node_pt][2 * ilayer - 1];
1369
1370 new_el_pt[1][ilayer]->node_pt(j_local + 9) =
1371 connected_node_pt[new_node_pt][2 * ilayer];
1372 new_el_pt[0][ilayer]->node_pt(j_local + 9) =
1373 connected_node_pt[new_node_pt][2 * ilayer];
1374
1375 new_el_pt[1][ilayer]->node_pt(j_local + 18) =
1376 connected_node_pt[new_node_pt][2 * ilayer + 1];
1377 new_el_pt[0][ilayer]->node_pt(j_local + 18) =
1378 connected_node_pt[new_node_pt][2 * ilayer + 1];
1379 }
1380
1381 // Nodes 6 and 7 in all elements are the same as nodes 2 and 5
1382 //------------------------------------------------------------
1383 // in previous element around the patch
1384 //-------------------------------------
1385 for (unsigned ilayer = 0; ilayer < nlayer; ilayer++)
1386 {
1387 for (unsigned j = 0; j < 3; j++)
1388 {
1389 unsigned offset = 9 * j;
1390 new_el_pt[2][ilayer]->node_pt(6 + offset) =
1391 new_el_pt[1][ilayer]->node_pt(2 + offset);
1392
1393 new_el_pt[1][ilayer]->node_pt(6 + offset) =
1394 new_el_pt[0][ilayer]->node_pt(2 + offset);
1395
1396 new_el_pt[0][ilayer]->node_pt(6 + offset) =
1397 new_el_pt[2][ilayer]->node_pt(2 + offset);
1398
1399 new_el_pt[2][ilayer]->node_pt(7 + offset) =
1400 new_el_pt[1][ilayer]->node_pt(5 + offset);
1401
1402 new_el_pt[1][ilayer]->node_pt(7 + offset) =
1403 new_el_pt[0][ilayer]->node_pt(5 + offset);
1404
1405 new_el_pt[0][ilayer]->node_pt(7 + offset) =
1406 new_el_pt[2][ilayer]->node_pt(5 + offset);
1407 }
1408 }
1409
1410 // Outer boundary is the last one
1411 Outer_boundary_id = maxb - 1;
1412
1413 // Number of identified in/outflow domain boundaries
1414 // (remember they're broken up into separeate boundaries with oomph-lib)
1415 unsigned nb_in_out = is_on_in_out_boundary.size();
1416 In_out_boundary_id.resize(nb_in_out);
1417
1418 // Now loop over the elements in the stacks again
1419 // and add all connected nodes to the appopriate non-FSI
1420 // boundary
1421 for (unsigned j_stack = 0; j_stack < 3; j_stack++)
1422 {
1423 // Bottom element
1424 FiniteElement* el_pt = new_el_pt[j_stack][0];
1425
1426 // Loop over nodes in bottom layer
1427 for (unsigned j = 0; j < 9; j++)
1428 {
1429 // Get nodes above...
1430 Node* nod_pt = el_pt->node_pt(j);
1431 Vector<Node*> layer_node_pt = connected_node_pt[nod_pt];
1432 unsigned n = layer_node_pt.size();
1433
1434 // Get boundary affiliation
1435 std::set<unsigned>* bnd_pt;
1436 nod_pt->get_boundaries_pt(bnd_pt);
1437
1438 // Loop over boundaries
1439 for (std::set<unsigned>::iterator it = (*bnd_pt).begin();
1440 it != (*bnd_pt).end();
1441 it++)
1442 {
1443 // Ignore FSI surface!
1444 if (!is_on_fsi_boundary[(*it)])
1445 {
1446 // Loop over connnected nodes in layers above
1447 unsigned ilayer = 0;
1448 for (unsigned k = 0; k < n; k++)
1449 {
1450 // Add to boundary
1451 add_boundary_node((*it), layer_node_pt[k]);
1452
1453 int face_index = 0;
1454
1455 // Use edge node on bottom node layer to assess
1456 // the element/s affiliation with boundary
1457 if (j == 1) face_index = -2;
1458 if (j == 3) face_index = -1;
1459 if (j == 5) face_index = 1;
1460 if (j == 7) face_index = 2;
1461
1462 if (face_index != 0)
1463 {
1464 // Use middle level in vertical direction
1465 // to assess the element's affiliation with boundary
1466 if (k % 2 == 1)
1467 {
1468 Boundary_element_pt[(*it)].push_back(
1469 new_el_pt[j_stack][ilayer]);
1470 Face_index_at_boundary[(*it)].push_back(face_index);
1471 ilayer++;
1472 }
1473 }
1474
1475 // Add to lookup scheme that allows the nodes
1476 // associated with an identified macroscopic in/outflow
1477 // boundary to recovered collectively.
1478
1479 // Loop over macroscopic in/outflow boundaries
1480 for (unsigned jj = 0; jj < nb_in_out; jj++)
1481 {
1482 if (is_on_in_out_boundary[jj][(*it)])
1483 {
1484 in_out_boundary_id_set[jj].insert((*it));
1485 }
1486 }
1487 }
1488 }
1489 }
1490
1491 // Last connected node is on outer boundary
1492 add_boundary_node(Outer_boundary_id, layer_node_pt[n - 1]);
1493
1494 // Use central node on bottom node layer to assess
1495 // the element/s affiliation with outer boundary
1496 if (j == 4)
1497 {
1498 Boundary_element_pt[Outer_boundary_id].push_back(
1499 new_el_pt[j_stack][nlayer - 1]);
1500 int face_index = 3;
1501 Face_index_at_boundary[Outer_boundary_id].push_back(face_index);
1502 }
1503 }
1504 }
1505
1506 // Cleanup
1507 delete face_el_pt;
1508 }
1509 }
1510
1511 // Copy boundary IDs across
1512 for (unsigned jj = 0; jj < n; jj++)
1513 {
1514 for (std::set<unsigned>::iterator it = in_out_boundary_id_set[jj].begin();
1515 it != in_out_boundary_id_set[jj].end();
1516 it++)
1517 {
1518 In_out_boundary_id[jj].push_back((*it));
1519 }
1520 }
1521
1522#ifdef PARANOID
1523 // Check
1524 unsigned nel = Element_pt.size();
1525 for (unsigned e = 0; e < nel; e++)
1526 {
1527 FiniteElement* el_pt = finite_element_pt(e);
1528 for (unsigned j = 0; j < 27; j++)
1529 {
1530 if (el_pt->node_pt(j) == 0)
1531 {
1532 // Throw an error
1533 std::ostringstream error_stream;
1534 error_stream << "Null node in element " << e << " node " << j
1535 << std::endl;
1536 throw OomphLibError(error_stream.str(),
1537 OOMPH_CURRENT_FUNCTION,
1538 OOMPH_EXCEPTION_LOCATION);
1539 }
1540 }
1541 }
1542#endif
1543
1544 // Average unit normals
1545 std::ofstream outfile;
1546 bool doc_normals = false; // keep alive for future debugging
1547 if (doc_normals) outfile.open("normals.dat");
1548 for (std::map<Node*, Vector<Vector<double>>>::iterator it = normals.begin();
1549 it != normals.end();
1550 it++)
1551 {
1552 Vector<double> unit_normal(3, 0.0);
1553 unsigned nnormals = ((*it).second).size();
1554 for (unsigned j = 0; j < nnormals; j++)
1555 {
1556 for (unsigned i = 0; i < 3; i++)
1557 {
1558 unit_normal[i] += ((*it).second)[j][i];
1559 }
1560 }
1561 double norm = 0.0;
1562 for (unsigned i = 0; i < 3; i++)
1563 {
1564 norm += unit_normal[i] * unit_normal[i];
1565 }
1566 for (unsigned i = 0; i < 3; i++)
1567 {
1568 unit_normal[i] /= sqrt(norm);
1569 }
1570
1571 Node* base_node_pt = (*it).first;
1572 Vector<double> base_pos(3);
1573 base_node_pt->position(base_pos);
1574 double h_thick;
1575 Thickness_fct_pt(base_pos, h_thick);
1576 Vector<Node*> layer_node_pt = connected_node_pt[base_node_pt];
1577 unsigned n = layer_node_pt.size();
1578 for (unsigned j = 0; j < n; j++)
1579 {
1580 for (unsigned i = 0; i < 3; i++)
1581 {
1582 layer_node_pt[j]->x(i) =
1583 base_pos[i] + h_thick * double(j + 1) / double(n) * unit_normal[i];
1584 }
1585 }
1586 if (doc_normals)
1587 {
1588 outfile << ((*it).first)->x(0) << " " << ((*it).first)->x(1) << " "
1589 << ((*it).first)->x(2) << " " << unit_normal[0] << " "
1590 << unit_normal[1] << " " << unit_normal[2] << "\n";
1591 }
1592 }
1593 if (doc_normals) outfile.close();
1594
1595 // Lookup scheme has now been setup yet
1596 Lookup_for_elements_next_boundary_is_setup = true;
1597 }
1598
1599} // namespace oomph
1600
1601#endif
e
e
Definition cfortran.h:571
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
oomph::Edge
Edge class.
Definition mesh.h:2700
oomph::FaceElement
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition elements.h:4342
oomph::FiniteElement
A general Finite Element class.
Definition elements.h:1317
oomph::FiniteElement::position
void position(const Vector< double > &zeta, Vector< double > &r) const
Return the parametrised position of the FiniteElement in its incarnation as a GeomObject,...
Definition elements.h:2680
oomph::FiniteElement::interpolated_zeta
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition elements.cc:4705
oomph::FiniteElement::size
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition elements.cc:4320
oomph::FiniteElement::construct_node
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition elements.h:2513
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition elements.cc:3992
oomph::FiniteElement::construct_boundary_node
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition elements.h:2542
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
oomph::Mesh
A general mesh class.
Definition mesh.h:67
oomph::Mesh::add_boundary_node
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
oomph::Mesh::Node_pt
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition mesh.h:183
oomph::Mesh::Boundary_element_pt
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,...
Definition mesh.h:91
oomph::Mesh::Lookup_for_elements_next_boundary_is_setup
bool Lookup_for_elements_next_boundary_is_setup
Flag to indicate that the lookup schemes for elements that are adjacent to the boundaries has been se...
Definition mesh.h:87
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition mesh.h:477
oomph::Mesh::Face_index_at_boundary
Vector< Vector< int > > Face_index_at_boundary
For the e-th finite element on boundary b, this is the index of the face that lies along that boundar...
Definition mesh.h:95
oomph::Mesh::set_nboundary
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition mesh.h:509
oomph::Mesh::set_boundary_coordinate_exists
void set_boundary_coordinate_exists(const unsigned &i)
Set boundary coordinate on the i-th boundary to be existing.
Definition mesh.h:584
oomph::Mesh::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition mesh.h:186
oomph::Node
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition nodes.h:906
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition oomph_definitions.h:229
oomph::SolidFiniteElement
SolidFiniteElement class.
Definition elements.h:3565
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
oomph::ThinLayerBrickOnTetMesh::Thickness_fct_pt
ThicknessFctPt Thickness_fct_pt
Function pointer to function that specifies the wall thickness as a fct of the coordinates of the inn...
Definition thin_layer_brick_on_tet_mesh.h:110
oomph::ThinLayerBrickOnTetMesh::Outer_boundary_id
unsigned Outer_boundary_id
Boundary ID of the "outer" surface – the non-wetted tube surface at a distance h_thick from the FSI s...
Definition thin_layer_brick_on_tet_mesh.h:102
oomph::ThinLayerBrickOnTetMesh::In_out_boundary_id
Vector< Vector< unsigned > > In_out_boundary_id
Vector of vectors containing the ids of the oomph-lib mesh boundaries that make up the specified in/o...
Definition thin_layer_brick_on_tet_mesh.h:106
oomph::ThinLayerBrickOnTetMesh::in_out_boundary_id
Vector< unsigned > in_out_boundary_id(const unsigned &boundary_id)
Access function to the vector containing the ids of the oomph-lib mesh boundaries that make up the sp...
Definition thin_layer_brick_on_tet_mesh.h:89
oomph::ThinLayerBrickOnTetMesh::FSI_boundary_id
Vector< unsigned > FSI_boundary_id
Vector of oomph-lib boundary ids that make up boundary on which the mesh was erected (typically the F...
Definition thin_layer_brick_on_tet_mesh.h:98
oomph::ThinLayerBrickOnTetMesh::ThinLayerBrickOnTetMesh
ThinLayerBrickOnTetMesh(Mesh *tet_mesh_pt, const Vector< unsigned > &boundary_ids, ThicknessFctPt thickness_fct_pt, const unsigned &nlayer, const Vector< Vector< unsigned > > &in_out_boundary_id, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Specify (quadratic) tet mesh, boundary IDs of boundary on which the current mesh is to b...
Definition thin_layer_brick_on_tet_mesh.template.cc:54
oomph::TimeStepper
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition timesteppers.h:231
oomph::Vector
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition Vector.h:58
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30