Public Types | Public Member Functions | Private Attributes | List of all members
oomph::ThinLayerBrickOnTetMesh< ELEMENT > Class Template Reference

#include <thin_layer_brick_on_tet_mesh.template.h>

+ Inheritance diagram for oomph::ThinLayerBrickOnTetMesh< ELEMENT >:

Public Types

typedef void(* ThicknessFctPt) (const Vector< double > &x, double &h_thick)
 
- Public Types inherited from oomph::Mesh
typedef void(FiniteElement::* SteadyExactSolutionFctPt) (const Vector< double > &x, Vector< double > &soln)
 
typedef void(FiniteElement::* UnsteadyExactSolutionFctPt) (const double &time, const Vector< double > &x, Vector< double > &soln)
 

Public Member Functions

 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)
 
Vector< unsigned > fsi_boundary_id ()
 
unsigned outer_boundary_id ()
 
Vector< unsigned > in_out_boundary_id (const unsigned &boundary_id)
 
- Public Member Functions inherited from oomph::BrickMeshBase
 BrickMeshBase ()
 Constructor (empty) More...
 
 BrickMeshBase (const BrickMeshBase &)=delete
 Broken copy constructor. More...
 
virtual ~BrickMeshBase ()
 Broken assignment operator. More...
 
void setup_boundary_element_info ()
 
void setup_boundary_element_info (std::ostream &outfile)
 
- Public Member Functions inherited from oomph::Mesh
 Mesh ()
 Default constructor. More...
 
 Mesh (const Vector< Mesh * > &sub_mesh_pt)
 
void merge_meshes (const Vector< Mesh * > &sub_mesh_pt)
 
virtual void reset_boundary_element_info (Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned >> &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
 Virtual function to perform the reset boundary elements info rutines. More...
 
template<class BULK_ELEMENT >
void doc_boundary_coordinates (const unsigned &b, std::ofstream &the_file)
 
virtual void scale_mesh (const double &factor)
 
 Mesh (const Mesh &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const Mesh &)=delete
 Broken assignment operator. More...
 
virtual ~Mesh ()
 Virtual Destructor to clean up all memory. More...
 
void flush_element_and_node_storage ()
 
void flush_element_storage ()
 
void flush_node_storage ()
 
Node *& node_pt (const unsigned long &n)
 Return pointer to global node n. More...
 
Node * node_pt (const unsigned long &n) const
 Return pointer to global node n (const version) More...
 
GeneralisedElement *& element_pt (const unsigned long &e)
 Return pointer to element e. More...
 
GeneralisedElement * element_pt (const unsigned long &e) const
 Return pointer to element e (const version) More...
 
const Vector< GeneralisedElement * > & element_pt () const
 Return reference to the Vector of elements. More...
 
Vector< GeneralisedElement * > & element_pt ()
 Return reference to the Vector of elements. More...
 
FiniteElement * finite_element_pt (const unsigned &e) const
 
Node *& boundary_node_pt (const unsigned &b, const unsigned &n)
 Return pointer to node n on boundary b. More...
 
Node * boundary_node_pt (const unsigned &b, const unsigned &n) const
 Return pointer to node n on boundary b. More...
 
void set_nboundary (const unsigned &nbound)
 Set the number of boundaries in the mesh. More...
 
void remove_boundary_nodes ()
 Clear all pointers to boundary nodes. More...
 
void remove_boundary_nodes (const unsigned &b)
 
void remove_boundary_node (const unsigned &b, Node *const &node_pt)
 Remove a node from the boundary b. More...
 
void add_boundary_node (const unsigned &b, Node *const &node_pt)
 Add a (pointer to) a node to the b-th boundary. More...
 
void copy_boundary_node_data_from_nodes ()
 
bool boundary_coordinate_exists (const unsigned &i) const
 Indicate whether the i-th boundary has an intrinsic coordinate. More...
 
unsigned long nelement () const
 Return number of elements in the mesh. More...
 
unsigned long nnode () const
 Return number of nodes in the mesh. More...
 
unsigned ndof_types () const
 Return number of dof types in mesh. More...
 
unsigned elemental_dimension () const
 Return number of elemental dimension in mesh. More...
 
unsigned nodal_dimension () const
 Return number of nodal dimension in mesh. More...
 
void add_node_pt (Node *const &node_pt)
 Add a (pointer to a) node to the mesh. More...
 
void add_element_pt (GeneralisedElement *const &element_pt)
 Add a (pointer to) an element to the mesh. More...
 
virtual void node_update (const bool &update_all_solid_nodes=false)
 
virtual void reorder_nodes (const bool &use_old_ordering=true)
 
virtual void get_node_reordering (Vector< Node * > &reordering, const bool &use_old_ordering=true) const
 
template<class BULK_ELEMENT , template< class > class FACE_ELEMENT>
void build_face_mesh (const unsigned &b, Mesh *const &face_mesh_pt)
 
unsigned self_test ()
 Self-test: Check elements and nodes. Return 0 for OK. More...
 
void max_and_min_element_size (double &max_size, double &min_size)
 
double total_size ()
 
void check_inverted_elements (bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
 
void check_inverted_elements (bool &mesh_has_inverted_elements)
 
unsigned check_for_repeated_nodes (const double &epsilon=1.0e-12)
 
Vector< Node * > prune_dead_nodes ()
 
unsigned nboundary () const
 Return number of boundaries. More...
 
unsigned long nboundary_node (const unsigned &ibound) const
 Return number of nodes on a particular boundary. More...
 
FiniteElement * boundary_element_pt (const unsigned &b, const unsigned &e) const
 Return pointer to e-th finite element on boundary b. More...
 
Node * get_some_non_boundary_node () const
 
unsigned nboundary_element (const unsigned &b) const
 Return number of finite elements that are adjacent to boundary b. More...
 
int face_index_at_boundary (const unsigned &b, const unsigned &e) const
 
virtual void dump (std::ofstream &dump_file, const bool &use_old_ordering=true) const
 Dump the data in the mesh into a file for restart. More...
 
void dump (const std::string &dump_file_name, const bool &use_old_ordering=true) const
 Dump the data in the mesh into a file for restart. More...
 
virtual void read (std::ifstream &restart_file)
 Read solution from restart file. More...
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
void output (std::ostream &outfile)
 Output for all elements. More...
 
void output (std::ostream &outfile, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output (FILE *file_pt)
 Output for all elements (C-style output) More...
 
void output (FILE *file_pt, const unsigned &nplot)
 Output at f(n_plot) points in each element (C-style output) More...
 
void output (const std::string &output_filename)
 Output for all elements. More...
 
void output (const std::string &output_filename, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt)
 Output a given Vector function at f(n_plot) points in each element. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt)
 
void output_boundaries (std::ostream &outfile)
 Output the nodes on the boundaries (into separate tecplot zones) More...
 
void output_boundaries (const std::string &output_filename)
 
void assign_initial_values_impulsive ()
 Assign initial values for an impulsive start. More...
 
void shift_time_values ()
 
void calculate_predictions ()
 
void set_nodal_and_elemental_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void set_mesh_level_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void set_consistent_pinned_values_for_continuation (ContinuationStorageScheme *const &continuation_stepper_pt)
 Set consistent values for pinned data in continuation. More...
 
bool does_pointer_correspond_to_mesh_data (double *const &parameter_pt)
 Does the double pointer correspond to any mesh data. More...
 
void set_nodal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 Set the timestepper associated with the nodal data in the mesh. More...
 
void set_elemental_internal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void compute_norm (double &norm)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Returns the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
bool is_mesh_distributed () const
 Boolean to indicate if Mesh has been distributed. More...
 
OomphCommunicator * communicator_pt () const
 
void delete_all_external_storage ()
 Wipe the storage for all externally-based elements. More...
 

Private Attributes

Vector< unsigned > FSI_boundary_id
 
unsigned Outer_boundary_id
 
Vector< Vector< unsigned > > In_out_boundary_id
 
ThicknessFctPt Thickness_fct_pt
 

Additional Inherited Members

- Static Public Attributes inherited from oomph::Mesh
static Steady< 0 > Default_TimeStepper
 The Steady Timestepper. More...
 
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
 Static boolean flag to control warning about mesh level timesteppers. More...
 
- Protected Member Functions inherited from oomph::Mesh
unsigned long assign_global_eqn_numbers (Vector< double * > &Dof_pt)
 Assign (global) equation numbers to the nodes. More...
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign local equation numbers in all elements. More...
 
void convert_to_boundary_node (Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
 
void convert_to_boundary_node (Node *&node_pt)
 
- Protected Attributes inherited from oomph::Mesh
Vector< Vector< Node * > > Boundary_node_pt
 
bool Lookup_for_elements_next_boundary_is_setup
 
Vector< Vector< FiniteElement * > > Boundary_element_pt
 
Vector< Vector< int > > Face_index_at_boundary
 
Vector< Node * > Node_pt
 Vector of pointers to nodes. More...
 
Vector< GeneralisedElement * > Element_pt
 Vector of pointers to generalised elements. More...
 
std::vector< bool > Boundary_coordinate_exists
 

Detailed Description

template<class ELEMENT>
class oomph::ThinLayerBrickOnTetMesh< ELEMENT >

Brick mesh layer built on top of a given tet mesh. Typically used in FSI problems where the tet mesh is the fluid mesh and this mesh acts as the solid mesh that surrounds the FSI interface.

Member Typedef Documentation

◆ ThicknessFctPt

template<class ELEMENT >
typedef void(* oomph::ThinLayerBrickOnTetMesh< ELEMENT >::ThicknessFctPt) (const Vector< double > &x, double &h_thick)

Function pointer to function that specifies the wall thickness as a fct of the coordinates of the inner surface

Constructor & Destructor Documentation

◆ ThinLayerBrickOnTetMesh()

template<class ELEMENT >
oomph::ThinLayerBrickOnTetMesh< ELEMENT >::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 be erected (in an FSI context this boundary tends to be the FSI boundary of the fluid mesh. Also specify the uniform thickness of layer, and the number of element layers. The vectors stored in in_out_boundary_ids contain the boundary IDs of the other boundaries in the tet mesh. In an FSI context these typically identify the in/outflow boundaries in the fluid mesh. The boundary enumeration of the current mesh follows the one of the underlying fluid mesh: The enumeration of the FSI boundary matches (to enable the setup of the FSI matching); the "in/outflow" faces in this mesh inherit the same enumeration as the in/outflow faces in the underlying fluid mesh. Finally, the "outer" boundary gets its own boundary ID. Timestepper defaults to steady pseudo-timestepper.

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  }
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::Mesh::add_boundary_node
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
oomph::Mesh::Node_pt
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
oomph::Mesh::Boundary_element_pt
Vector< Vector< FiniteElement * > > Boundary_element_pt
Definition: mesh.h:91
oomph::Mesh::Lookup_for_elements_next_boundary_is_setup
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:87
oomph::Mesh::Boundary_coordinate_exists
std::vector< bool > Boundary_coordinate_exists
Definition: mesh.h:190
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
oomph::Mesh::Face_index_at_boundary
Vector< Vector< int > > Face_index_at_boundary
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::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
oomph::ThinLayerBrickOnTetMesh::in_out_boundary_id
Vector< unsigned > in_out_boundary_id(const unsigned &boundary_id)
Definition: thin_layer_brick_on_tet_mesh.template.h:90
oomph::ThinLayerBrickOnTetMesh::Thickness_fct_pt
ThicknessFctPt Thickness_fct_pt
Definition: thin_layer_brick_on_tet_mesh.template.h:112
oomph::ThinLayerBrickOnTetMesh::Outer_boundary_id
unsigned Outer_boundary_id
Definition: thin_layer_brick_on_tet_mesh.template.h:104
oomph::ThinLayerBrickOnTetMesh::In_out_boundary_id
Vector< Vector< unsigned > > In_out_boundary_id
Definition: thin_layer_brick_on_tet_mesh.template.h:108
oomph::ThinLayerBrickOnTetMesh::FSI_boundary_id
Vector< unsigned > FSI_boundary_id
Definition: thin_layer_brick_on_tet_mesh.template.h:100
Eigen::Vector
Matrix< Type, Size, 1 > Vector
\cpp11 Size×1 vector of type Type.
Definition: Eigen/Eigen/src/Core/Matrix.h:515
s
RealScalar s
Definition: level1_cplx_impl.h:130
nb
int nb
Definition: level2_impl.h:286
k
char char char int int * k
Definition: level2_impl.h:374
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
Global_Parameters::is_on_fsi_boundary
bool is_on_fsi_boundary(Node *nod_pt)
Boolean to identify if node is on fsi boundary.
Definition: unstructured_two_d_fsi.cc:70
plotDoE.x
list x
Definition: plotDoE.py:28
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::Mesh::add_boundary_node(), b, oomph::Mesh::Boundary_coordinate_exists, oomph::Mesh::Boundary_element_pt, oomph::Mesh::boundary_element_pt(), oomph::Mesh::boundary_node_pt(), e(), oomph::Mesh::Element_pt, oomph::Mesh::element_pt(), oomph::Mesh::Face_index_at_boundary, oomph::Mesh::face_index_at_boundary(), oomph::Mesh::finite_element_pt(), oomph::ThinLayerBrickOnTetMesh< ELEMENT >::FSI_boundary_id, oomph::Node::get_boundaries_pt(), i, oomph::ThinLayerBrickOnTetMesh< ELEMENT >::in_out_boundary_id(), oomph::ThinLayerBrickOnTetMesh< ELEMENT >::In_out_boundary_id, oomph::FaceElement::interpolated_x(), oomph::FiniteElement::interpolated_zeta(), Global_Parameters::is_on_fsi_boundary(), j, k, oomph::Mesh::Lookup_for_elements_next_boundary_is_setup, n, nb, oomph::Mesh::nboundary_element(), oomph::Mesh::nboundary_node(), oomph::Edge::node1_pt(), oomph::Edge::node2_pt(), oomph::FiniteElement::node_pt(), oomph::Mesh::Node_pt, oomph::FaceElement::normal_sign(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::ThinLayerBrickOnTetMesh< ELEMENT >::Outer_boundary_id, oomph::FaceElement::outer_unit_normal(), oomph::Node::position(), s, oomph::FaceElement::set_boundary_number_in_bulk_mesh(), oomph::Node::set_coordinates_on_boundary(), oomph::Mesh::set_nboundary(), size, sqrt(), oomph::ThinLayerBrickOnTetMesh< ELEMENT >::Thickness_fct_pt, plotDoE::x, oomph::Node::x(), and Eigen::zeta().

Member Function Documentation

◆ fsi_boundary_id()

template<class ELEMENT >
Vector<unsigned> oomph::ThinLayerBrickOnTetMesh< ELEMENT >::fsi_boundary_id ( )
inline

Access functions to the Vector of oomph-lib boundary ids that make up boundary on which the mesh was erected (typically the FSI interface in an FSI problem)

75  {
76  return FSI_boundary_id;
77  }

References oomph::ThinLayerBrickOnTetMesh< ELEMENT >::FSI_boundary_id.

◆ in_out_boundary_id()

template<class ELEMENT >
Vector<unsigned> oomph::ThinLayerBrickOnTetMesh< ELEMENT >::in_out_boundary_id ( const unsigned &  boundary_id)
inline

Access function to the vector containing the ids of the oomph-lib mesh boundaries that make up the specified in/outflow boundaries as specified in constructor.

91  {
92  return In_out_boundary_id[boundary_id];
93  }

References oomph::ThinLayerBrickOnTetMesh< ELEMENT >::In_out_boundary_id.

Referenced by oomph::ThinLayerBrickOnTetMesh< ELEMENT >::ThinLayerBrickOnTetMesh().

◆ outer_boundary_id()

template<class ELEMENT >
unsigned oomph::ThinLayerBrickOnTetMesh< ELEMENT >::outer_boundary_id ( )
inline

Boundary ID of the "outer" surface – in an FSI context this is the non-wetted tube surface at a distance h_thick from the FSI surface

83  {
84  return Outer_boundary_id;
85  }

References oomph::ThinLayerBrickOnTetMesh< ELEMENT >::Outer_boundary_id.

Member Data Documentation

◆ FSI_boundary_id

template<class ELEMENT >
Vector<unsigned> oomph::ThinLayerBrickOnTetMesh< ELEMENT >::FSI_boundary_id
private

Vector of oomph-lib boundary ids that make up boundary on which the mesh was erected (typically the FSI interface in an FSI problem)

Referenced by oomph::ThinLayerBrickOnTetMesh< ELEMENT >::fsi_boundary_id(), and oomph::ThinLayerBrickOnTetMesh< ELEMENT >::ThinLayerBrickOnTetMesh().

◆ In_out_boundary_id

template<class ELEMENT >
Vector<Vector<unsigned> > oomph::ThinLayerBrickOnTetMesh< ELEMENT >::In_out_boundary_id
private

Vector of vectors containing the ids of the oomph-lib mesh boundaries that make up the specified in/outflow boundaries

Referenced by oomph::ThinLayerBrickOnTetMesh< ELEMENT >::in_out_boundary_id(), and oomph::ThinLayerBrickOnTetMesh< ELEMENT >::ThinLayerBrickOnTetMesh().

◆ Outer_boundary_id

template<class ELEMENT >
unsigned oomph::ThinLayerBrickOnTetMesh< ELEMENT >::Outer_boundary_id
private

Boundary ID of the "outer" surface – the non-wetted tube surface at a distance h_thick from the FSI surface

Referenced by oomph::ThinLayerBrickOnTetMesh< ELEMENT >::outer_boundary_id(), and oomph::ThinLayerBrickOnTetMesh< ELEMENT >::ThinLayerBrickOnTetMesh().

◆ Thickness_fct_pt

template<class ELEMENT >
ThicknessFctPt oomph::ThinLayerBrickOnTetMesh< ELEMENT >::Thickness_fct_pt
private

Function pointer to function that specifies the wall thickness as a fct of the coordinates of the inner surface

Referenced by oomph::ThinLayerBrickOnTetMesh< ELEMENT >::ThinLayerBrickOnTetMesh().

The documentation for this class was generated from the following files: