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