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