tetgen_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_TETGEN_MESH_TEMPLATE_CC
31 #define OOMPH_TETGEN_MESH_TEMPLATE_CC
32 
33 
34 #include<algorithm>
35 
36 #include "tetgen_mesh.template.h"
37 #include "../generic/Telements.h"
38 #include "../generic/map_matrix.h"
39 
40 
41 
42 namespace oomph
43 {
44 
45 
46 ///////////////////////////////////////////////////////////////////////
47 ///////////////////////////////////////////////////////////////////////
48 ///////////////////////////////////////////////////////////////////////
49 
50 
51 
52 //========================================================================
53 /// Build unstructured tet mesh based on output from scaffold
54 //========================================================================
55 template <class ELEMENT>
57  const bool &use_attributes)
58 {
59  // Mesh can only be built with 3D Telements.
60  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
61 
62  // Create space for elements
63  unsigned nelem=Tmp_mesh_pt->nelement();
64  Element_pt.resize(nelem);
65 
66  // Create space for nodes
67  unsigned nnode_scaffold=Tmp_mesh_pt->nnode();
68  Node_pt.resize(nnode_scaffold);
69 
70  // Set number of boundaries
71  unsigned nbound=Tmp_mesh_pt->nboundary();
72  set_nboundary(nbound);
73  Boundary_element_pt.resize(nbound);
74  Face_index_at_boundary.resize(nbound);
75 
76  //If we have different regions, then resize the region
77  //information
78  if(use_attributes)
79  {
80  Boundary_region_element_pt.resize(nbound);
81  Face_index_region_at_boundary.resize(nbound);
82  }
83 
84  // Build elements
85  for (unsigned e=0;e<nelem;e++)
86  {
87  Element_pt[e]=new ELEMENT;
88  }
89 
90  // Number of nodes per element
91  unsigned nnod_el=Tmp_mesh_pt->finite_element_pt(0)->nnode();
92 
93  // Setup map to check the (pseudo-)global node number
94  // Nodes whose number is zero haven't been copied across
95  // into the mesh yet.
96  std::map<Node*,unsigned> global_number;
97  unsigned global_count=0;
98 
99  // Map of element attribute pairs
100  std::map<double,Vector<FiniteElement*> > element_attribute_map;
101 
102  // Loop over elements in scaffold mesh, visit their nodes
103  for (unsigned e=0;e<nelem;e++)
104  {
105  // Loop over all nodes in element
106  for (unsigned j=0;j<nnod_el;j++)
107  {
108 
109  // Pointer to node in the scaffold mesh
110  Node* scaffold_node_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(j);
111 
112  // Get the (pseudo-)global node number in scaffold mesh
113  // (It's zero [=default] if not visited this one yet)
114  unsigned j_global=global_number[scaffold_node_pt];
115 
116  // Haven't done this one yet
117  if (j_global==0)
118  {
119  // Get pointer to set of mesh boundaries that this
120  // scaffold node occupies; NULL if the node is not on any boundary
121  std::set<unsigned>* boundaries_pt;
122  scaffold_node_pt->get_boundaries_pt(boundaries_pt);
123 
124  // Is it on boundaries?
125  if (boundaries_pt!=0)
126  {
127  // Create new boundary node
128  Node* new_node_pt=finite_element_pt(e)->
129  construct_boundary_node(j,time_stepper_pt);
130 
131  // Give it a number (not necessarily the global node
132  // number in the scaffold mesh -- we just need something
133  // to keep track...)
134  global_count++;
135  global_number[scaffold_node_pt]=global_count;
136 
137  // Add to boundaries
138  for(std::set<unsigned>::iterator it=boundaries_pt->begin();
139  it!=boundaries_pt->end();++it)
140  {
141  add_boundary_node(*it,new_node_pt);
142  }
143  }
144  // Build normal node
145  else
146  {
147  // Create new normal node
148  finite_element_pt(e)->construct_node(j,time_stepper_pt);
149 
150  // Give it a number (not necessarily the global node
151  // number in the scaffold mesh -- we just need something
152  // to keep track...)
153  global_count++;
154  global_number[scaffold_node_pt]=global_count;
155  }
156 
157  // Copy new node, created using the NEW element's construct_node
158  // function into global storage, using the same global
159  // node number that we've just associated with the
160  // corresponding node in the scaffold mesh
161  Node_pt[global_count-1]=finite_element_pt(e)->node_pt(j);
162 
163  // Assign coordinates
164  Node_pt[global_count-1]->x(0)=scaffold_node_pt->x(0);
165  Node_pt[global_count-1]->x(1)=scaffold_node_pt->x(1);
166  Node_pt[global_count-1]->x(2)=scaffold_node_pt->x(2);
167 
168  }
169  // This one has already been done: Copy across
170  else
171  {
172  finite_element_pt(e)->node_pt(j)=Node_pt[j_global-1];
173  }
174  }
175 
176  //Store the attributes in the map
177  if(use_attributes)
178  {
179  element_attribute_map[Tmp_mesh_pt->element_attribute(e)].push_back(
180  finite_element_pt(e));
181  }
182  }
183 
184  //Now let's construct lists
185  //Find the number of attributes
186  if(use_attributes)
187  {
188  unsigned n_attribute = element_attribute_map.size();
189  //There are n_attribute different regions
190  Region_element_pt.resize(n_attribute);
191  Region_attribute.resize(n_attribute);
192  //Copy the vectors in the map over to our internal storage
193  unsigned count = 0;
194  for(std::map<double,Vector<FiniteElement*> >::iterator it =
195  element_attribute_map.begin(); it != element_attribute_map.end();++it)
196  {
197  Region_attribute[count] = it->first;
198  Region_element_pt[count] = it->second;
199  ++count;
200  }
201  }
202 
203  // At this point we've created all the elements and
204  // created their vertex nodes. Now we need to create
205  // the additional (midside and internal) nodes!
206  unsigned boundary_id;
207 
208  // Get number of nodes along element edge and dimension of element (3)
209  // from first element
210  unsigned n_node_1d=finite_element_pt(0)->nnode_1d();
211 
212  // At the moment we're only able to deal with nnode_1d=2 or 3.
213  if ((n_node_1d!=2)&&(n_node_1d!=3))
214  {
215  std::ostringstream error_message;
216  error_message << "Mesh generation from tetgen currently only works\n";
217  error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
218  error_message << "for nnode_1d=" << n_node_1d << std::endl;
219 
220  throw OomphLibError(error_message.str(),
221  OOMPH_CURRENT_FUNCTION,
222  OOMPH_EXCEPTION_LOCATION);
223  }
224 
225  // Spatial dimension of element = number of local coordinates
226  unsigned dim=finite_element_pt(0)->dim();
227 
228  // Storage for the local coordinate of the new node
229  Vector<double> s(dim);
230 
231  // Get number of nodes in the element from first element
232  unsigned n_node=finite_element_pt(0)->nnode();
233 
234  // Storage for each global edge of the mesh
235  unsigned n_global_edge = Tmp_mesh_pt->nglobal_edge();
236  Vector<Vector<Node*> > nodes_on_global_edge(n_global_edge);
237 
238  //Storage for each global face of the mesh
239  unsigned n_global_face = Tmp_mesh_pt->nglobal_face();
240  Vector<Vector<Node*> > nodes_on_global_face(n_global_face);
241 
242  // Map storing the mid-side of an edge; edge identified by
243  // pointers to vertex nodes in scaffold mesh
244  //MapMatrix<Node*,Node*> central_edge_node_pt;
245  //Node* edge_node1_pt=0;
246  //Node* edge_node2_pt=0;
247 
248  // Map storing the mid point of a face; face identified by
249  // set of pointers to vertex nodes in scaffold mesh
250  //std::map<std::set<Node*>,Node*> central_face_node_pt;
251  //std::set<Node*> face_nodes_pt;
252 
253  //Mapping of Tetgen faces to face nodes in the enriched element
254  unsigned face_map[4] = {1,0,2,3};
255 
256  //Storage for the faces shared by the edges
257  const unsigned faces_on_edge[6][2]={
258  {0,1},{0,2},{1,2},{0,3},{2,3},{1,3}};
259 
260  // Loop over all elements
261  for(unsigned e=0;e<nelem;e++)
262  {
263  //Cache pointers to the elements
264  FiniteElement* const elem_pt = this->finite_element_pt(e);
265  FiniteElement* const tmp_elem_pt = Tmp_mesh_pt->finite_element_pt(e);
266 
267  // The number of edge nodes is 4 + 6*(n_node1d-2)
268  unsigned n_edge_node = 4 + 6*(n_node_1d-2);
269 
270  //Now loop over the edge nodes
271  //assuming that the numbering scheme is the same as the triangles
272  //which puts edge nodes in ascending order.
273  //We don't have any higher than quadratic at the moment, so it's
274  //a bit academic really.
275 
276  //Only bother if we actually have some edge nodes
277  if(n_edge_node > 4)
278  {
279  //Start from node number 4
280  unsigned n=4;
281 
282  //Loop over the edges
283  for(unsigned j=0;j<6;++j)
284  {
285  //Find the global edge index
286  unsigned edge_index = Tmp_mesh_pt->edge_index(e,j);
287 
288  //Use the intersection of the appropriate faces to determine
289  //whether the boundaries on which an edge lies
290  std::set<unsigned> edge_boundaries;
291  for(unsigned i=0;i<2;++i)
292  {
293  unsigned face_boundary_id =
294  Tmp_mesh_pt->face_boundary(e,faces_on_edge[j][i]);
295  if(face_boundary_id > 0)
296  {
297  edge_boundaries.insert(face_boundary_id);
298  }
299  }
300 
301  //If the nodes on the edge have not been allocated, construct them
302  if(nodes_on_global_edge[edge_index].size()==0)
303  {
304  //Now loop over the nodes on the edge
305  for(unsigned j2=0;j2<n_node_1d-2;++j2)
306  {
307  //Storage for the new node
308  Node* new_node_pt = 0;
309 
310  //If the edge is on a boundary, determined from the
311  //scaffold mesh, construct a boundary node
312  //The use of the scaffold mesh's edge_boundary data structure
313  //ensures that a boundary node is created, even if the faces of
314  //the current element do not lie on boundaries
315  if(Tmp_mesh_pt->edge_boundary(edge_index) == true)
316  {
317  new_node_pt =
318  elem_pt->construct_boundary_node(n,time_stepper_pt);
319  //Add it to the boundaries in the set,
320  //remembering to subtract one to get to the oomph-lib numbering
321  //scheme
322  for(std::set<unsigned>::iterator it = edge_boundaries.begin();
323  it!=edge_boundaries.end();++it)
324  {
325  this->add_boundary_node((*it)-1,new_node_pt);
326  }
327  }
328  //Otherwise construct a normal node
329  else
330  {
331  new_node_pt =
332  elem_pt->construct_node(n,time_stepper_pt);
333  }
334 
335  //Find the local coordinates of the node
336  elem_pt->local_coordinate_of_node(n,s);
337 
338  //Find the coordinates of the new node from the exiting and
339  //fully-functional element in the scaffold mesh
340  for(unsigned i=0;i<dim;++i)
341  {
342  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s,i);
343  }
344 
345  //Add the newly created node to the global node list
346  Node_pt.push_back(new_node_pt);
347  //Add to the edge index
348  nodes_on_global_edge[edge_index].push_back(new_node_pt);
349  //Increment the local node number
350  ++n;
351  } //end of loop over edge nodes
352  }
353  //Otherwise just set the pointers (orientation assumed the same)
354  else
355  {
356  for(unsigned j2=0;j2<n_node_1d-2;++j2)
357  {
358  elem_pt->node_pt(n) =
359  nodes_on_global_edge[edge_index][j2];
360  //It is possible that the edge may be on additional boundaries
361  //through another element
362  //So add again (note that this function will not add to
363  //boundaries twice)
364  for(std::set<unsigned>::iterator it = edge_boundaries.begin();
365  it!=edge_boundaries.end();++it)
366  {
367  this->add_boundary_node((*it)-1,elem_pt->node_pt(n));
368  }
369  ++n;
370  }
371  }
372  } //End of loop over edges
373 
374  //Deal with enriched elements
375  if(n_node == 15)
376  {
377  //Now loop over the faces
378  for(unsigned j=0;j<4;++j)
379  {
380  //Find the boundary id of the face (need to map between node numbers
381  //and the face)
382  boundary_id = Tmp_mesh_pt->face_boundary(e,face_map[j]);
383 
384  //Find the global face index (need to map between node numbers and
385  //the face)
386  unsigned face_index = Tmp_mesh_pt->face_index(e,face_map[j]);
387 
388  //If the nodes on the face have not been allocated
389  if(nodes_on_global_face[face_index].size()==0)
390  {
391  //Storage for the new node
392  Node* new_node_pt=0;
393 
394  //If the face is on a boundary, construct a boundary node
395  if(boundary_id > 0)
396  {
397  new_node_pt =
398  elem_pt->construct_boundary_node(n,time_stepper_pt);
399  //Add it to the boundary
400  this->add_boundary_node(boundary_id-1,new_node_pt);
401  }
402  //Otherwise construct a normal node
403  else
404  {
405  new_node_pt =
406  elem_pt->construct_node(n,time_stepper_pt);
407  }
408 
409  //Find the local coordinates of the node
410  elem_pt->local_coordinate_of_node(n,s);
411 
412  //Find the coordinates of the new node from the exiting and
413  //fully-functional element in the scaffold mesh
414  for(unsigned i=0;i<dim;++i)
415  {
416  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s,i);
417  }
418 
419  //Add the newly created node to the global node list
420  Node_pt.push_back(new_node_pt);
421  //Add to the face index
422  nodes_on_global_face[face_index].push_back(new_node_pt);
423  //Increment the local node number
424  ++n;
425  }
426  //Otherwise just set the single node from the face element
427  else
428  {
429  elem_pt->node_pt(n) = nodes_on_global_face[face_index][0];
430  ++n;
431  }
432  } //end of loop over faces
433 
434  //Construct the element's central node, which is not on a boundary
435  {
436  Node* new_node_pt=
437  finite_element_pt(e)->construct_node(n,time_stepper_pt);
438  Node_pt.push_back(new_node_pt);
439 
440  // Find the local coordinates of the node
441  elem_pt->local_coordinate_of_node(n,s);
442 
443  // Find the coordinates of the new node from the existing
444  // and fully-functional element in the scaffold mesh
445  for(unsigned i=0;i<dim;i++)
446  {
447  new_node_pt->x(i)=tmp_elem_pt->interpolated_x(s,i);
448  }
449  }
450  } //End of enriched case
451 
452  } //End of case for edge nodes
453 
454 
455  //Now loop over the faces to setup the information about elements
456  //adjacent to the boundary
457  for(unsigned j=0;j<4;++j)
458  {
459  //Find the boundary id of the face
460  boundary_id = Tmp_mesh_pt->face_boundary(e,j);
461 
462  if(boundary_id > 0)
463  {
464  Boundary_element_pt[boundary_id-1].push_back(elem_pt);
465  //Need to put a shift in here because of an inconsistent naming
466  //convention between tetgen and our faces
467  //Tetgen Face 0 is our Face 3
468  //Tetgen Face 1 is our Face 2
469  //Tetgen Face 2 is our Face 1
470  //Tetgen Face 3 is our Face 0
471  Face_index_at_boundary[boundary_id-1].push_back(3-j);
472 
473  //If using regions set up the boundary information
474  if(use_attributes)
475  {
476  //Element adjacent to boundary
477  Boundary_region_element_pt[boundary_id-1]
478  [static_cast<unsigned>(Tmp_mesh_pt->element_attribute(e))].
479  push_back(elem_pt);
480  //Need to put a shift in here because of an inconsistent naming
481  //convention between triangle and face elements
482  Face_index_region_at_boundary[boundary_id-1]
483  [static_cast<unsigned>(Tmp_mesh_pt->element_attribute(e))].
484  push_back(3-j);
485  }
486  }
487  } //End of loop over faces
488 
489 
490  //Lookup scheme has now been setup
491  Lookup_for_elements_next_boundary_is_setup=true;
492 
493 
494  /*
495 
496  //For standard quadratic elements all nodes are edge nodes
497  unsigned n_vertex_and_edge_node = n_node;
498  //If we have enriched elements, there are only 10 vertex and edge nodes
499  if(n_node==15)
500  {
501  //There are only 10 vertex and edge nodes
502  n_vertex_and_edge_node = 10;
503  }
504 
505  // Loop over the new (edge) nodes in the element and create them.
506  for(unsigned j=4;j<n_vertex_and_edge_node;j++)
507  {
508 
509  // Figure out edge nodes
510  switch(j)
511  {
512 
513  // Node 4 is between nodes 0 and 1
514  case 4:
515 
516  edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
517  edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
518  break;
519 
520 
521  // Node 5 is between nodes 0 and 2
522  case 5:
523 
524  edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
525  edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
526  break;
527 
528  // Node 6 is between nodes 0 and 3
529  case 6:
530 
531  edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
532  edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
533  break;
534 
535  // Node 7 is between nodes 1 and 2
536  case 7:
537 
538  edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
539  edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
540  break;
541 
542  // Node 8 is between nodes 2 and 3
543  case 8:
544 
545  edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
546  edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
547  break;
548 
549  // Node 9 is between nodes 1 and 3
550  case 9:
551 
552  edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
553  edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
554  break;
555 
556  break;
557 
558  //Error
559  default:
560 
561  throw OomphLibError("More than ten edge nodes in Tet Element",
562  OOMPH_CURRENT_FUNCTION,
563  OOMPH_EXCEPTION_LOCATION);
564  }
565 
566 
567 
568 
569  // Do we need a boundary node?
570  bool need_boundary_node=false;
571 
572  // hierher At some point fine tune via intersection of boundary sets
573  if (edge_node1_pt->is_on_boundary()&&edge_node2_pt->is_on_boundary())
574  {
575  need_boundary_node=true;
576  }
577 
578  // Do we need a new node?
579  if (central_edge_node_pt(edge_node1_pt,edge_node2_pt)==0)
580  {
581  Node* new_node_pt=0;
582 
583  // Create new boundary node
584  if (need_boundary_node)
585  {
586  new_node_pt=finite_element_pt(e)->
587  construct_boundary_node(j,time_stepper_pt);
588  }
589  // Create new normal node
590  else
591  {
592  new_node_pt=finite_element_pt(e)->
593  construct_node(j,time_stepper_pt);
594  }
595  Node_pt.push_back(new_node_pt);
596 
597  // Now indicate existence of newly created mideside node in map
598  central_edge_node_pt(edge_node1_pt,edge_node2_pt)=new_node_pt;
599  central_edge_node_pt(edge_node2_pt,edge_node1_pt)=new_node_pt;
600 
601  // What are the node's local coordinates?
602  finite_element_pt(e)->local_coordinate_of_node(j,s);
603 
604  // Find the coordinates of the new node from the existing
605  // and fully-functional element in the scaffold mesh
606  for(unsigned i=0;i<dim;i++)
607  {
608  new_node_pt->x(i)=
609  Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
610  }
611  }
612  else
613  {
614  // Set pointer to the existing node
615  finite_element_pt(e)->node_pt(j)=
616  central_edge_node_pt(edge_node1_pt,edge_node2_pt);
617  }
618 
619  } // end of loop over new nodes
620 
621  //Need to sort out the face nodes
622  if(n_node==15)
623  {
624  // Loop over the new (face) nodes in the element and create them.
625  for(unsigned j=10;j<14;j++)
626  {
627  //Empty the set of face nodes
628  face_nodes_pt.clear();
629  // Figure out face nodes
630  switch(j)
631  {
632 
633  // Node 10 is between nodes 0 and 1 and 3
634  case 10:
635 
636  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
637  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
638  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
639  break;
640 
641  // Node 11 is between nodes 0 and 1 and 2
642  case 11:
643 
644  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
645  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
646  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
647  break;
648 
649  // Node 12 is between nodes 0 and 2 and 3
650  case 12:
651 
652  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
653  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
654  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
655  break;
656 
657  // Node 13 is between nodes 1 and 2 and 3
658  case 13:
659 
660  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
661  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
662  face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
663  break;
664 
665 
666  //Error
667  default:
668 
669  throw OomphLibError("More than four face nodes in Tet Element",
670  OOMPH_CURRENT_FUNCTION,
671  OOMPH_EXCEPTION_LOCATION);
672  }
673 
674  // Do we need a boundary node?
675  bool need_boundary_node=false;
676 
677  //Work it out from the face boundary
678  boundary_id = Tmp_mesh_pt->face_boundary(e,face_map[j-10]);
679  //If it's non-zero then we do need to create a boundary node
680  if(boundary_id!=0) {need_boundary_node=true;}
681 
682  // Do we need a new node?
683  if (central_face_node_pt[face_nodes_pt]==0)
684  {
685  Node* new_node_pt=0;
686 
687  // Create a new boundary node
688  if (need_boundary_node)
689  {
690  new_node_pt=finite_element_pt(e)->
691  construct_boundary_node(j,time_stepper_pt);
692  //Add it to the boundary
693  add_boundary_node(boundary_id-1,new_node_pt);
694  }
695  // Create new normal node
696  else
697  {
698  new_node_pt=finite_element_pt(e)->
699  construct_node(j,time_stepper_pt);
700  }
701  Node_pt.push_back(new_node_pt);
702 
703  // Now indicate existence of newly created mideside node in map
704  central_face_node_pt[face_nodes_pt]=new_node_pt;
705 
706  // What are the node's local coordinates?
707  finite_element_pt(e)->local_coordinate_of_node(j,s);
708 
709  // Find the coordinates of the new node from the existing
710  // and fully-functional element in the scaffold mesh
711  for(unsigned i=0;i<dim;i++)
712  {
713  new_node_pt->x(i)=
714  Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
715  }
716  }
717  else
718  {
719  // Set pointer to the existing node
720  finite_element_pt(e)->node_pt(j)=
721  central_face_node_pt[face_nodes_pt];
722  }
723  } //End of loop over face nodes
724 
725  //Construct the element's central node, which is not on a boundary
726  {
727  Node* new_node_pt=
728  finite_element_pt(e)->construct_node(14,time_stepper_pt);
729  Node_pt.push_back(new_node_pt);
730 
731  // What are the node's local coordinates?
732  finite_element_pt(e)->local_coordinate_of_node(14,s);
733  // Find the coordinates of the new node from the existing
734  // and fully-functional element in the scaffold mesh
735  for(unsigned i=0;i<dim;i++)
736  {
737  new_node_pt->x(i)=
738  Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
739  }
740  }
741  } //End of enriched case
742 
743  } //end of loop over elements
744 
745 
746  //Boundary conditions
747 
748  // Matrix map to check if a node has already been added to
749  // the boundary number b
750  MapMatrixMixed<Node*,unsigned,bool> bound_node_done;
751 
752  // Loop over elements
753  for (unsigned e=0;e<nelem;e++)
754  {
755  // Loop over new local nodes
756  for(unsigned j=4;j<n_node;j++)
757  {
758  // Pointer to the element's local node
759  Node* loc_node_pt=finite_element_pt(e)->node_pt(j);
760 
761  // This will have to be changed for higher-order elements
762  //=======================================================
763 
764  // These are the face nodes on the element's face 0:
765  if ( (j==4) || (j==5) || (j==7) )
766  {
767  boundary_id=Tmp_mesh_pt->face_boundary(e,0);
768  if (boundary_id!=0)
769  {
770  if (!bound_node_done(loc_node_pt,boundary_id-1))
771  {
772  add_boundary_node(boundary_id-1,loc_node_pt);
773  bound_node_done(loc_node_pt,boundary_id-1)=true;
774  }
775  }
776  }
777 
778 
779  // These are the face nodes on the element's face 1:
780  if ( (j==4) || (j==6) || (j==9) )
781  {
782  boundary_id=Tmp_mesh_pt->face_boundary(e,1);
783  if (boundary_id!=0)
784  {
785  if (!bound_node_done(loc_node_pt,boundary_id-1))
786  {
787  add_boundary_node(boundary_id-1,loc_node_pt);
788  bound_node_done(loc_node_pt,boundary_id-1)=true;
789  }
790  }
791  }
792 
793  // These are the face nodes on the element's face 2:
794  if ( (j==5) || (j==6) || (j==8) )
795  {
796  boundary_id=Tmp_mesh_pt->face_boundary(e,2);
797  if (boundary_id!=0)
798  {
799  if (!bound_node_done(loc_node_pt,boundary_id-1))
800  {
801  add_boundary_node(boundary_id-1,loc_node_pt);
802  bound_node_done(loc_node_pt,boundary_id-1)=true;
803  }
804  }
805  }
806 
807 
808  // These are the face nodes on the element's face 3:
809  if ( (j==7) || (j==8) || (j==9) )
810  {
811  boundary_id=Tmp_mesh_pt->face_boundary(e,3);
812  if (boundary_id!=0)
813  {
814  if (!bound_node_done(loc_node_pt,boundary_id-1))
815  {
816  add_boundary_node(boundary_id-1,loc_node_pt);
817  bound_node_done(loc_node_pt,boundary_id-1)=true;
818  }
819  }
820  }
821 
822  } //end for j
823 
824  */
825 
826  } //end for e
827 
828 
829 } // end function
830 
831 
832 //=========================================================================
833 /// Transfer tetgenio data from the input to the output
834 /// The output is assumed to have been constructed and "empty"
835  template<class ELEMENT>
836  void TetgenMesh<ELEMENT>::deep_copy_of_tetgenio(tetgenio* const &input_pt,
837  tetgenio* &output_pt)
838  {
839  //Copy data values
840  output_pt->firstnumber = input_pt->firstnumber;
841  output_pt->mesh_dim = input_pt->mesh_dim;
842  output_pt->useindex = input_pt->useindex;
843 
844  //Copy the number of points
845  output_pt->numberofpoints = input_pt->numberofpoints;
846  output_pt->numberofpointattributes = input_pt->numberofpointattributes;
847  output_pt->numberofpointmtrs = input_pt->numberofpointmtrs;
848 
849  const int n_point = output_pt->numberofpoints;
850  if(n_point > 0)
851  {
852  output_pt->pointlist = new double[n_point * 3];
853  //Copy the values
854  for(int n=0;n<(n_point * 3);++n)
855  {
856  output_pt->pointlist[n] = input_pt->pointlist[n];
857  }
858 
859  //If there are point markers copy as well
860  if(input_pt->pointmarkerlist != (int*) NULL)
861  {
862  output_pt->pointmarkerlist = new int[n_point];
863  for(int n=0;n<n_point;++n)
864  {
865  output_pt->pointmarkerlist[n] = input_pt->pointmarkerlist[n];
866  }
867  }
868 
869  const int n_attr = output_pt->numberofpointattributes;
870  if(n_attr > 0)
871  {
872  output_pt->pointattributelist = new double[n_point * n_attr];
873  for(int n=0;n<(n_point * n_attr);++n)
874  {
875  output_pt->pointattributelist[n] =
876  input_pt->pointattributelist[n];
877  }
878  }
879  } //End of point case
880 
881  //Now add in metric tensors (if there are any)
882  const int n_point_mtr = output_pt->numberofpointmtrs;
883  if(n_point_mtr > 0)
884  {
885  output_pt->pointmtrlist = new double[n_point * n_point_mtr];
886  for(int n=0;n<(n_point * n_point_mtr);++n)
887  {
888  output_pt->pointmtrlist[n] = input_pt->pointmtrlist[n];
889  }
890  }
891 
892  //Next tetrahedrons
893  output_pt->numberoftetrahedra = input_pt->numberoftetrahedra;
894  output_pt->numberofcorners = input_pt->numberofcorners;
895  output_pt->numberoftetrahedronattributes =
896  input_pt->numberoftetrahedronattributes;
897 
898  const int n_tetra = output_pt->numberoftetrahedra;
899  const int n_corner = output_pt->numberofcorners;
900  if(n_tetra > 0)
901  {
902  output_pt->tetrahedronlist = new int[n_tetra * n_corner];
903  for(int n=0;n<(n_tetra * n_corner);++n)
904  {
905  output_pt->tetrahedronlist[n] = input_pt->tetrahedronlist[n];
906  }
907 
908  //Add in the volume constriants
909  if(input_pt->tetrahedronvolumelist != (double*) NULL)
910  {
911  output_pt->tetrahedronvolumelist = new double[n_tetra];
912  for(int n=0;n<n_tetra;++n)
913  {
914  output_pt->tetrahedronvolumelist[n] =
915  input_pt->tetrahedronvolumelist[n];
916  }
917  }
918 
919  //Add in the attributes
920  const int n_tetra_attr = output_pt->numberoftetrahedronattributes;
921  if(n_tetra_attr > 0)
922  {
923  output_pt->tetrahedronattributelist = new double[n_tetra * n_tetra_attr];
924  for(int n=0;n<(n_tetra * n_tetra_attr);++n)
925  {
926  output_pt->tetrahedronattributelist[n] =
927  input_pt->tetrahedronattributelist[n];
928  }
929  }
930 
931  //Add in the neighbour list
932  if(input_pt->neighborlist != (int*) NULL)
933  {
934  output_pt->neighborlist = new int[n_tetra * 4];
935  for(int n=0;n<(n_tetra * 4);++n)
936  {
937  output_pt->neighborlist = input_pt->neighborlist;
938  }
939  }
940  } //End of tetrahedron section
941 
942  //Now arrange the facet list
943  output_pt->numberoffacets = input_pt->numberoffacets;
944  const int n_facet = output_pt->numberoffacets;
945  if(n_facet > 0)
946  {
947  output_pt->facetlist = new tetgenio::facet[n_facet];
948  for(int n=0;n<n_facet;++n)
949  {
950  tetgenio::facet *input_f_pt = &input_pt->facetlist[n];
951  tetgenio::facet *output_f_pt = &output_pt->facetlist[n];
952 
953  //Copy polygons and holes from the facets
954  output_f_pt->numberofpolygons = input_f_pt->numberofpolygons;
955 
956  //Loop over polygons
957  const int n_poly = output_f_pt->numberofpolygons;
958  if(n_poly > 0)
959  {
960  output_f_pt->polygonlist = new tetgenio::polygon[n_poly];
961  for(int p=0;p<n_poly;++p)
962  {
963  tetgenio::polygon *output_p_pt = &output_f_pt->polygonlist[p];
964  tetgenio::polygon *input_p_pt = &input_f_pt->polygonlist[p];
965  //Find out how many vertices each polygon has
966  output_p_pt->numberofvertices = input_p_pt->numberofvertices;
967  //Now copy of the vertices
968  const int n_vertex = output_p_pt->numberofvertices;
969  if(n_vertex > 0)
970  {
971  output_p_pt->vertexlist = new int[n_vertex];
972  for(int v=0;v<n_vertex;++v)
973  {
974  output_p_pt->vertexlist[v] = input_p_pt->vertexlist[v];
975  }
976  }
977  }
978  }
979 
980  //Hole information
981  output_f_pt->numberofholes = input_f_pt->numberofholes;
982  const int n_hole = output_f_pt->numberofholes;
983  if(n_hole > 0)
984  {
985  throw OomphLibError("Don't know how to handle holes yet\n",
986  OOMPH_CURRENT_FUNCTION,
987  OOMPH_EXCEPTION_LOCATION);
988  }
989  } //end of loop over facets
990 
991  //Add the facetmarkers if there are any
992  if(input_pt->facetmarkerlist != (int*) NULL)
993  {
994  output_pt->facetmarkerlist = new int[n_facet];
995  for(int n=0;n<n_facet;++n)
996  {
997  output_pt->facetmarkerlist[n] = input_pt->facetmarkerlist[n];
998  }
999  }
1000  }
1001 
1002  //Now the holes
1003  output_pt->numberofholes = input_pt->numberofholes;
1004  const int n_hole = output_pt->numberofholes;
1005  if(n_hole > 0)
1006  {
1007  output_pt->holelist = new double[n_hole * 3];
1008  for(int n=0;n<(n_hole * 3);++n)
1009  {
1010  output_pt->holelist[n] = input_pt->holelist[n];
1011  }
1012  }
1013 
1014  //Now the regions
1015  output_pt->numberofregions = input_pt->numberofregions;
1016  const int n_region = output_pt->numberofregions;
1017  if(n_region > 0)
1018  {
1019  output_pt->regionlist = new double[n_region * 5];
1020  for(int n=0;n<(n_region * 5);++n)
1021  {
1022  output_pt->regionlist[n] = input_pt->regionlist[n];
1023  }
1024  }
1025 
1026  //Now the facet constraints
1027  output_pt->numberoffacetconstraints = input_pt->numberoffacetconstraints;
1028  const int n_facet_const = output_pt->numberoffacetconstraints;
1029  if(n_facet_const > 0)
1030  {
1031  output_pt->facetconstraintlist = new double[n_facet_const * 2];
1032  for(int n=0;n<(n_facet_const * 2);++n)
1033  {
1034  output_pt->facetconstraintlist[n] = input_pt->facetconstraintlist[n];
1035  }
1036  }
1037 
1038  //Now the segment constraints
1039  output_pt->numberofsegmentconstraints = input_pt->numberofsegmentconstraints;
1040  const int n_seg_const = output_pt->numberofsegmentconstraints;
1041  if(n_seg_const > 0)
1042  {
1043  output_pt->segmentconstraintlist = new double[n_seg_const * 2];
1044  for(int n=0;n<(n_seg_const * 2);++n)
1045  {
1046  output_pt->segmentconstraintlist[n] = input_pt->segmentconstraintlist[n];
1047  }
1048  }
1049 
1050  //Now the face list
1051  output_pt->numberoftrifaces = input_pt->numberoftrifaces;
1052  const int n_tri_face = output_pt->numberoftrifaces;
1053  if(n_tri_face > 0)
1054  {
1055  output_pt->trifacelist = new int[n_tri_face * 3];
1056  for(int n=0;n<(n_tri_face * 3);++n)
1057  {
1058  output_pt->trifacelist[n] = input_pt->trifacelist[n];
1059  }
1060 
1061  output_pt->trifacemarkerlist = new int[n_tri_face];
1062  for(int n=0;n<n_tri_face;++n)
1063  {
1064  output_pt->trifacemarkerlist[n] = input_pt->trifacemarkerlist[n];
1065  }
1066  }
1067 
1068  //Now the edge list
1069  output_pt->numberofedges = input_pt->numberofedges;
1070  const int n_edge = output_pt->numberofedges;
1071  if(n_edge > 0)
1072  {
1073  output_pt->edgelist = new int[n_edge * 2];
1074  for(int n=0;n<(n_edge * 2);++n)
1075  {
1076  output_pt->edgelist[n] = input_pt->edgelist[n];
1077  }
1078 
1079  output_pt->edgemarkerlist = new int[n_edge];
1080  for(int n=0;n<n_edge;++n)
1081  {
1082  output_pt->edgemarkerlist[n] = input_pt->edgemarkerlist[n];
1083  }
1084  }
1085  }
1086 
1087 
1088 
1089 //======================================================================
1090 /// Setup boundary coordinate on boundary b which is
1091 /// assumed to be planar. Boundary coordinates are the
1092 /// x-y coordinates in the plane of that boundary with the
1093 /// x-axis along the line from the (lexicographically)
1094 /// "lower left" to the "upper right" node. The y axis
1095 /// is obtained by taking the cross-product of the positive
1096 /// x direction with the outer unit normal computed by
1097 /// the face elements (or its negative if switch_normal is set
1098 /// to true). Doc faces in output file.
1099 //======================================================================
1100 template <class ELEMENT>
1102  const unsigned& b,
1103  const bool& switch_normal,
1104  std::ofstream& outfile)
1105 {
1106 
1107  // Temporary storage for face elements
1108  Vector<FiniteElement*> face_el_pt;
1109 
1110  // Loop over all elements on boundaries
1111  unsigned nel=this->nboundary_element(b);
1112  if (nel>0)
1113  {
1114  // Loop over the bulk elements adjacent to boundary b
1115  for(unsigned e=0;e<nel;e++)
1116  {
1117  // Get pointer to the bulk element that is adjacent to boundary b
1118  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b,e);
1119 
1120  //Find the index of the face of element e along boundary b
1121  int face_index = this->face_index_at_boundary(b,e);
1122 
1123  // Create new face element
1124  face_el_pt.push_back(new DummyFaceElement<ELEMENT>(
1125  bulk_elem_pt,face_index));
1126 
1127  // Output faces?
1128  if (outfile.is_open())
1129  {
1130  face_el_pt[face_el_pt.size()-1]->output(outfile);
1131  }
1132  }
1133 
1134 
1135  // Loop over all nodes to find the lower left and upper right ones
1136  Node* lower_left_node_pt=this->boundary_node_pt(b,0);
1137  Node* upper_right_node_pt=this->boundary_node_pt(b,0);
1138 
1139  // Loop over all nodes on the boundary
1140  unsigned nnod=this->nboundary_node(b);
1141  for (unsigned j=0;j<nnod;j++)
1142  {
1143 
1144  //Get node
1145  Node* nod_pt=this->boundary_node_pt(b,j);
1146 
1147  // Primary criterion for lower left: Does it have a lower z-coordinate?
1148  if (nod_pt->x(2)<lower_left_node_pt->x(2))
1149  {
1150  lower_left_node_pt=nod_pt;
1151  }
1152  // ...or is it a draw?
1153  else if (nod_pt->x(2)==lower_left_node_pt->x(2))
1154  {
1155  // If it's a draw: Does it have a lower y-coordinate?
1156  if (nod_pt->x(1)<lower_left_node_pt->x(1))
1157  {
1158  lower_left_node_pt=nod_pt;
1159  }
1160  // ...or is it a draw?
1161  else if (nod_pt->x(1)==lower_left_node_pt->x(1))
1162  {
1163 
1164  // If it's a draw: Does it have a lower x-coordinate?
1165  if (nod_pt->x(0)<lower_left_node_pt->x(0))
1166  {
1167  lower_left_node_pt=nod_pt;
1168  }
1169  }
1170  }
1171 
1172  // Primary criterion for upper right: Does it have a higher z-coordinate?
1173  if (nod_pt->x(2)>upper_right_node_pt->x(2))
1174  {
1175  upper_right_node_pt=nod_pt;
1176  }
1177  // ...or is it a draw?
1178  else if (nod_pt->x(2)==upper_right_node_pt->x(2))
1179  {
1180  // If it's a draw: Does it have a higher y-coordinate?
1181  if (nod_pt->x(1)>upper_right_node_pt->x(1))
1182  {
1183  upper_right_node_pt=nod_pt;
1184  }
1185  // ...or is it a draw?
1186  else if (nod_pt->x(1)==upper_right_node_pt->x(1))
1187  {
1188 
1189  // If it's a draw: Does it have a higher x-coordinate?
1190  if (nod_pt->x(0)>upper_right_node_pt->x(0))
1191  {
1192  upper_right_node_pt=nod_pt;
1193  }
1194  }
1195  }
1196  }
1197 
1198  // Prepare storage for boundary coordinates
1199  Vector<double> zeta(2);
1200 
1201  // Unit vector connecting lower left and upper right nodes
1202  Vector<double> b0(3);
1203  b0[0]=upper_right_node_pt->x(0)-lower_left_node_pt->x(0);
1204  b0[1]=upper_right_node_pt->x(1)-lower_left_node_pt->x(1);
1205  b0[2]=upper_right_node_pt->x(2)-lower_left_node_pt->x(2);
1206 
1207  // Normalise
1208  double inv_length=1.0/sqrt(b0[0]*b0[0]+b0[1]*b0[1]+b0[2]*b0[2]);
1209  b0[0]*=inv_length;
1210  b0[1]*=inv_length;
1211  b0[2]*=inv_length;
1212 
1213  // Get (outer) unit normal to first face element
1214  Vector<double> normal(3);
1215  Vector<double> s(2,0.0);
1216  dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[0])->
1217  outer_unit_normal(s,normal);
1218 
1219  if (switch_normal)
1220  {
1221  normal[0]=-normal[0];
1222  normal[1]=-normal[1];
1223  normal[2]=-normal[2];
1224  }
1225 
1226 #ifdef PARANOID
1227 
1228  // Check that all elements have the same normal
1229  for(unsigned e=0;e<nel;e++)
1230  {
1231  // Get (outer) unit normal to face element
1232  Vector<double> my_normal(3);
1233  dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[0])->
1234  outer_unit_normal(s,my_normal);
1235 
1236  // Dot product should be one!
1237  double error=
1238  normal[0]*my_normal[0]+
1239  normal[1]*my_normal[1]+
1240  normal[2]*my_normal[2];
1241 
1242  error-=1.0;
1243  if (switch_normal) error+=1.0;
1244 
1245  if (error>Tolerance_for_boundary_finding)
1246  {
1247  std::ostringstream error_message;
1248  error_message
1249  << "Error in alingment of normals (dot product-1) "
1250  << error << " for element " << e << std::endl
1251  << "exceeds tolerance specified by the static member data\n"
1252  << "TetMeshBase::Tolerance_for_boundary_finding = "
1253  << Tolerance_for_boundary_finding << std::endl
1254  << "This usually means that the boundary is not planar.\n\n"
1255  << "You can suppress this error message by recompiling \n"
1256  << "recompiling without PARANOID or by changing the tolerance.\n";
1257  throw OomphLibError(error_message.str(),
1258  OOMPH_CURRENT_FUNCTION,
1259  OOMPH_EXCEPTION_LOCATION);
1260  }
1261  }
1262 #endif
1263 
1264 
1265 
1266 
1267 
1268  // Cross-product to get second in-plane vector, normal to b0
1269  Vector<double> b1(3);
1270  b1[0]=b0[1]*normal[2]-b0[2]*normal[1];
1271  b1[1]=b0[2]*normal[0]-b0[0]*normal[2];
1272  b1[2]=b0[0]*normal[1]-b0[1]*normal[0];
1273 
1274 
1275 
1276  // Assign boundary coordinates: projection onto the axes
1277  for (unsigned j=0;j<nnod;j++)
1278  {
1279  //Get node
1280  Node* nod_pt=this->boundary_node_pt(b,j);
1281 
1282  // Difference vector to lower left corner
1283  Vector<double> delta(3);
1284  delta[0]=nod_pt->x(0)-lower_left_node_pt->x(0);
1285  delta[1]=nod_pt->x(1)-lower_left_node_pt->x(1);
1286  delta[2]=nod_pt->x(2)-lower_left_node_pt->x(2);
1287 
1288  // Project
1289  zeta[0]=delta[0]*b0[0]+delta[1]*b0[1]+delta[2]*b0[2];
1290  zeta[1]=delta[0]*b1[0]+delta[1]*b1[1]+delta[2]*b1[2];
1291 
1292 #ifdef PARANOID
1293 
1294  // Check:
1295  double error=
1296  pow(lower_left_node_pt->x(0)+zeta[0]*b0[0]+zeta[1]*b1[0]-nod_pt->x(0),2)+
1297  pow(lower_left_node_pt->x(1)+zeta[0]*b0[1]+zeta[1]*b1[1]-nod_pt->x(1),2)+
1298  pow(lower_left_node_pt->x(2)+zeta[0]*b0[2]+zeta[1]*b1[2]-nod_pt->x(2),2);
1299 
1300  if (sqrt(error)>Tolerance_for_boundary_finding)
1301  {
1302  std::ostringstream error_message;
1303  error_message
1304  << "Error in setup of boundary coordinate "
1305  << sqrt(error) << std::endl
1306  << "exceeds tolerance specified by the static member data\n"
1307  << "TetMeshBase::Tolerance_for_boundary_finding = "
1308  << Tolerance_for_boundary_finding << std::endl
1309  << "This usually means that the boundary is not planar.\n\n"
1310  << "You can suppress this error message by recompiling \n"
1311  << "recompiling without PARANOID or by changing the tolerance.\n";
1312  throw OomphLibError(error_message.str(),
1313  OOMPH_CURRENT_FUNCTION,
1314  OOMPH_EXCEPTION_LOCATION);
1315  }
1316 #endif
1317 
1318  // Set boundary coordinate
1319  nod_pt->set_coordinates_on_boundary(b,zeta);
1320  }
1321 
1322  }
1323 
1324  // Indicate that boundary coordinate has been set up
1325  Boundary_coordinate_exists[b]=true;
1326 
1327  // Cleanup
1328  unsigned n=face_el_pt.size();
1329  for (unsigned e=0;e<n;e++)
1330  {
1331  delete face_el_pt[e];
1332  }
1333 
1334 }
1335 
1336 
1337 
1338 
1339 
1340 
1341 
1342 
1343 //======================================================================
1344 /// Snap boundaries specified by the IDs listed in boundary_id to
1345 /// a quadratric surface, specified in the file
1346 /// quadratic_surface_file_name. This is usually used with vmtk-based
1347 /// meshes for which oomph-lib's xda to poly conversion code produces the files
1348 /// "quadratic_fsi_boundary.dat" and "quadratic_outer_solid_boundary.dat"
1349 /// which specify the quadratic FSI boundary (for the fluid and the solid)
1350 /// and the quadratic representation of the outer boundary of the solid.
1351 /// When used with these files, the flag switch_normal should be
1352 /// set to true when calling the function for the outer boundary of the
1353 /// solid. The DocInfo object can be used to label optional output
1354 /// files. (Uses directory and label).
1355 //======================================================================
1356 template <class ELEMENT>
1358  const Vector<unsigned>& boundary_id,
1359  const std::string& quadratic_surface_file_name,
1360  const bool& switch_normal,
1361  DocInfo& doc_info)
1362 {
1363 
1364  // Aux storage for processing input
1365  char dummy[101];
1366 
1367  // Prepare to doc nodes that couldn't be snapped
1368  std::ofstream the_file_non_snapped;
1369  std::string non_snapped_filename="non_snapped_nodes_"+doc_info.label()+".dat";
1370 
1371  // Read the number of nodes and elements (quadratic facets)
1372  std::ifstream infile(quadratic_surface_file_name.c_str(),std::ios_base::in);
1373  unsigned n_node;
1374  infile >> n_node;
1375 
1376  // Ignore rest of line
1377  infile.getline(dummy, 100);
1378 
1379  // Number of quadratic facets
1380  unsigned nel;
1381  infile>> nel;
1382 
1383  // Ignore rest of line
1384  infile.getline(dummy, 100);
1385 
1386  // Check that the number of elements matches
1387  if (nel!=boundary_id.size())
1388  {
1389  std::ostringstream error_message;
1390  error_message
1391  << "Number of quadratic facets specified in "
1392  << quadratic_surface_file_name << "is: " << nel
1393  << "\nThis doesn't match the number of planar boundaries \n"
1394  << "specified in boundary_id which is " << boundary_id.size()
1395  << std::endl;
1396  throw OomphLibError(error_message.str(),
1397  OOMPH_CURRENT_FUNCTION,
1398  OOMPH_EXCEPTION_LOCATION);
1399  }
1400 
1401  // Temporary storage for face elements
1403 
1404  // Loop over quadratic face elements -- one for each facet
1405  for(unsigned e=0;e<nel;e++)
1406  {
1407  face_el_pt.push_back(new FreeStandingFaceElement<ELEMENT>);
1408  }
1409 
1410 
1411  // Now build nodes
1412  unsigned n_dim=3;
1413  unsigned n_position_type=1;
1414  unsigned initial_n_value=0;
1415  Vector<Node*> nod_pt(n_node);
1416  unsigned node_nmbr;
1417  std::map<unsigned,unsigned> node_number;
1418  std::ofstream node_file;
1419  for (unsigned j=0;j<n_node;j++)
1420  {
1421  nod_pt[j]=new BoundaryNode<Node>(n_dim,n_position_type,initial_n_value);
1422  infile >> nod_pt[j]->x(0);
1423  infile >> nod_pt[j]->x(1);
1424  infile >> nod_pt[j]->x(2);
1425  infile >> node_nmbr;
1426  node_number[node_nmbr]=j;
1427  }
1428 
1429 
1430  // Now assign nodes to elements -- each element represents
1431  // distinct boundary; assign enumeration as specified by
1432  // boundary_id.
1433  for(unsigned e=0;e<nel;e++)
1434  {
1435  unsigned nnod_el=face_el_pt[e]->nnode();
1436  unsigned j_global;
1437  for (unsigned j=0;j<nnod_el;j++)
1438  {
1439  infile >> j_global;
1440  face_el_pt[e]->node_pt(j)=nod_pt[node_number[j_global]];
1441  face_el_pt[e]->node_pt(j)->add_to_boundary(boundary_id[e]);
1442  }
1443  face_el_pt[e]->set_boundary_number_in_bulk_mesh(boundary_id[e]);
1444  face_el_pt[e]->set_nodal_dimension(3);
1445  }
1446 
1447 
1448  // Setup boundary coordinates for each facet, using
1449  // the same strategy as for the simplex boundaries
1450  // (there's some code duplication here but it doesn't
1451  // seem worth it to rationlise this as the interface would
1452  // be pretty clunky).
1453  for (unsigned e=0;e<nel;e++)
1454  {
1455  Vector<Vector<double> >vertex_boundary_coord(3);
1456 
1457  // Loop over all nodes to find the lower left and upper right ones
1458  Node* lower_left_node_pt=face_el_pt[e]->node_pt(0);
1459  Node* upper_right_node_pt=face_el_pt[e]->node_pt(0);
1460 
1461  // Loop over all vertex nodes
1462  for (unsigned j=0;j<3;j++)
1463  {
1464  //Get node
1465  Node* nod_pt=face_el_pt[e]->node_pt(j);
1466 
1467  // Primary criterion for lower left: Does it have a lower z-coordinate?
1468  if (nod_pt->x(2)<lower_left_node_pt->x(2))
1469  {
1470  lower_left_node_pt=nod_pt;
1471  }
1472  // ...or is it a draw?
1473  else if (nod_pt->x(2)==lower_left_node_pt->x(2))
1474  {
1475  // If it's a draw: Does it have a lower y-coordinate?
1476  if (nod_pt->x(1)<lower_left_node_pt->x(1))
1477  {
1478  lower_left_node_pt=nod_pt;
1479  }
1480  // ...or is it a draw?
1481  else if (nod_pt->x(1)==lower_left_node_pt->x(1))
1482  {
1483 
1484  // If it's a draw: Does it have a lower x-coordinate?
1485  if (nod_pt->x(0)<lower_left_node_pt->x(0))
1486  {
1487  lower_left_node_pt=nod_pt;
1488  }
1489  }
1490  }
1491 
1492  // Primary criterion for upper right: Does it have a higher z-coordinate?
1493  if (nod_pt->x(2)>upper_right_node_pt->x(2))
1494  {
1495  upper_right_node_pt=nod_pt;
1496  }
1497  // ...or is it a draw?
1498  else if (nod_pt->x(2)==upper_right_node_pt->x(2))
1499  {
1500  // If it's a draw: Does it have a higher y-coordinate?
1501  if (nod_pt->x(1)>upper_right_node_pt->x(1))
1502  {
1503  upper_right_node_pt=nod_pt;
1504  }
1505  // ...or is it a draw?
1506  else if (nod_pt->x(1)==upper_right_node_pt->x(1))
1507  {
1508 
1509  // If it's a draw: Does it have a higher x-coordinate?
1510  if (nod_pt->x(0)>upper_right_node_pt->x(0))
1511  {
1512  upper_right_node_pt=nod_pt;
1513  }
1514  }
1515  }
1516  }
1517 
1518  // Prepare storage for boundary coordinates
1519  Vector<double> zeta(2);
1520 
1521  // Unit vector connecting lower left and upper right nodes
1522  Vector<double> b0(3);
1523  b0[0]=upper_right_node_pt->x(0)-lower_left_node_pt->x(0);
1524  b0[1]=upper_right_node_pt->x(1)-lower_left_node_pt->x(1);
1525  b0[2]=upper_right_node_pt->x(2)-lower_left_node_pt->x(2);
1526 
1527  // Normalise
1528  double inv_length=1.0/sqrt(b0[0]*b0[0]+b0[1]*b0[1]+b0[2]*b0[2]);
1529  b0[0]*=inv_length;
1530  b0[1]*=inv_length;
1531  b0[2]*=inv_length;
1532 
1533  // Get (outer) unit normal to face element -- note that
1534  // with the enumeration chosen in oomph-lib's xda to poly
1535  // conversion code the sign below is correct for the outer
1536  // unit normal on the FSI interface.
1537  Vector<double> tang1(3);
1538  Vector<double> tang2(3);
1539  Vector<double> normal(3);
1540  tang1[0]=face_el_pt[e]->node_pt(1)->x(0)-face_el_pt[e]->node_pt(0)->x(0);
1541  tang1[1]=face_el_pt[e]->node_pt(1)->x(1)-face_el_pt[e]->node_pt(0)->x(1);
1542  tang1[2]=face_el_pt[e]->node_pt(1)->x(2)-face_el_pt[e]->node_pt(0)->x(2);
1543  tang2[0]=face_el_pt[e]->node_pt(2)->x(0)-face_el_pt[e]->node_pt(0)->x(0);
1544  tang2[1]=face_el_pt[e]->node_pt(2)->x(1)-face_el_pt[e]->node_pt(0)->x(1);
1545  tang2[2]=face_el_pt[e]->node_pt(2)->x(2)-face_el_pt[e]->node_pt(0)->x(2);
1546  normal[0]=-(tang1[1]*tang2[2]-tang1[2]*tang2[1]);
1547  normal[1]=-(tang1[2]*tang2[0]-tang1[0]*tang2[2]);
1548  normal[2]=-(tang1[0]*tang2[1]-tang1[1]*tang2[0]);
1549 
1550  // Normalise
1551  inv_length=
1552  1.0/sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
1553  normal[0]*=inv_length;
1554  normal[1]*=inv_length;
1555  normal[2]*=inv_length;
1556 
1557  // Change direction -- usually for outer boundary of solid
1558  if (switch_normal)
1559  {
1560  normal[0]=-normal[0];
1561  normal[1]=-normal[1];
1562  normal[2]=-normal[2];
1563  }
1564 
1565  // Cross-product to get second in-plane vector, normal to b0
1566  Vector<double> b1(3);
1567  b1[0]=b0[1]*normal[2]-b0[2]*normal[1];
1568  b1[1]=b0[2]*normal[0]-b0[0]*normal[2];
1569  b1[2]=b0[0]*normal[1]-b0[1]*normal[0];
1570 
1571  // Assign boundary coordinates for vertex nodes: projection onto the axes
1572  for (unsigned j=0;j<3;j++)
1573  {
1574  //Get node
1575  Node* nod_pt=face_el_pt[e]->node_pt(j);
1576 
1577  // Difference vector to lower left corner
1578  Vector<double> delta(3);
1579  delta[0]=nod_pt->x(0)-lower_left_node_pt->x(0);
1580  delta[1]=nod_pt->x(1)-lower_left_node_pt->x(1);
1581  delta[2]=nod_pt->x(2)-lower_left_node_pt->x(2);
1582 
1583  // Project
1584  zeta[0]=delta[0]*b0[0]+delta[1]*b0[1]+delta[2]*b0[2];
1585  zeta[1]=delta[0]*b1[0]+delta[1]*b1[1]+delta[2]*b1[2];
1586 
1587  vertex_boundary_coord[j].resize(2);
1588  vertex_boundary_coord[j][0]=zeta[0];
1589  vertex_boundary_coord[j][1]=zeta[1];
1590 
1591 #ifdef PARANOID
1592 
1593  // Check:
1594  double error=
1595  pow(lower_left_node_pt->x(0)+zeta[0]*b0[0]+zeta[1]*b1[0]-nod_pt->x(0),2)+
1596  pow(lower_left_node_pt->x(1)+zeta[0]*b0[1]+zeta[1]*b1[1]-nod_pt->x(1),2)+
1597  pow(lower_left_node_pt->x(2)+zeta[0]*b0[2]+zeta[1]*b1[2]-nod_pt->x(2),2);
1598 
1599  if (sqrt(error)>Tolerance_for_boundary_finding)
1600  {
1601  std::ostringstream error_message;
1602  error_message
1603  << "Error in setup of boundary coordinate "
1604  << sqrt(error) << std::endl
1605  << "exceeds tolerance specified by the static member data\n"
1606  << "TetMeshBase::Tolerance_for_boundary_finding = "
1607  << Tolerance_for_boundary_finding << std::endl
1608  << "This usually means that the boundary is not planar.\n\n"
1609  << "You can suppress this error message by recompiling \n"
1610  << "recompiling without PARANOID or by changing the tolerance.\n";
1611  throw OomphLibError(error_message.str(),
1612  OOMPH_CURRENT_FUNCTION,
1613  OOMPH_EXCEPTION_LOCATION);
1614  }
1615 #endif
1616 
1617  // Set boundary coordinate
1618  nod_pt->set_coordinates_on_boundary(boundary_id[e],zeta);
1619 
1620  }
1621 
1622  // Assign boundary coordinates of three midside nodes by linear
1623  // interpolation (corresponding to a flat facet)
1624 
1625  // Node 3 is between 0 and 1
1626  zeta[0]=0.5*(vertex_boundary_coord[0][0]+vertex_boundary_coord[1][0]);
1627  zeta[1]=0.5*(vertex_boundary_coord[0][1]+vertex_boundary_coord[1][1]);
1628  face_el_pt[e]->node_pt(3)->set_coordinates_on_boundary(boundary_id[e],zeta);
1629 
1630  // Node 4 is between 1 and 2
1631  zeta[0]=0.5*(vertex_boundary_coord[1][0]+vertex_boundary_coord[2][0]);
1632  zeta[1]=0.5*(vertex_boundary_coord[1][1]+vertex_boundary_coord[2][1]);
1633  face_el_pt[e]->node_pt(4)->set_coordinates_on_boundary(boundary_id[e],zeta);
1634 
1635  // Node 5 is between 2 and 0
1636  zeta[0]=0.5*(vertex_boundary_coord[2][0]+vertex_boundary_coord[0][0]);
1637  zeta[1]=0.5*(vertex_boundary_coord[2][1]+vertex_boundary_coord[0][1]);
1638  face_el_pt[e]->node_pt(5)->set_coordinates_on_boundary(boundary_id[e],zeta);
1639 
1640  }
1641 
1642 
1643  // Loop over elements/facets = boundaries to snap
1644  bool success=true;
1645  for(unsigned b=0;b<nel;b++)
1646  {
1647 
1648  //Doc boundary coordinates on quadratic patches
1649  std::ofstream the_file;
1650  std::ofstream the_file_before;
1651  std::ofstream the_file_after;
1652  if (doc_info.is_doc_enabled())
1653  {
1654  std::ostringstream filename;
1655  filename << doc_info.directory() << "/quadratic_coordinates_"
1656  << doc_info.label() << b << ".dat";
1657  the_file.open(filename.str().c_str());
1658 
1659  std::ostringstream filename1;
1660  filename1 << doc_info.directory() << "/quadratic_nodes_before_"
1661  << doc_info.label() << b << ".dat";
1662  the_file_before.open(filename1.str().c_str());
1663 
1664  std::ostringstream filename2;
1665  filename2 << doc_info.directory() << "/quadratic_nodes_after_"
1666  << doc_info.label() << b << ".dat";
1667  the_file_after.open(filename2.str().c_str());
1668  }
1669 
1670  //Cast the element pointer
1671  FreeStandingFaceElement<ELEMENT>* el_pt=face_el_pt[b];
1672 
1673  // Doc boundary coordinate on quadratic facet representation
1674  if (doc_info.is_doc_enabled())
1675  {
1676  Vector<double> s(2);
1677  Vector<double> zeta(2);
1678  Vector<double> x(3);
1679  unsigned n_plot=3;
1680  the_file << el_pt->tecplot_zone_string(n_plot);
1681 
1682  // Loop over plot points
1683  unsigned num_plot_points=el_pt->nplot_points(n_plot);
1684  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1685  {
1686  // Get local coordinates of plot point
1687  el_pt->get_s_plot(iplot,n_plot,s);
1688  el_pt->interpolated_zeta(s,zeta);
1689  el_pt->interpolated_x(s,x);
1690  for (unsigned i=0;i<3;i++)
1691  {
1692  the_file << x[i] << " ";
1693  }
1694  for (unsigned i=0;i<2;i++)
1695  {
1696  the_file << zeta[i] << " ";
1697  }
1698  the_file << std::endl;
1699  }
1700  el_pt->write_tecplot_zone_footer(the_file,n_plot);
1701 
1702 // std::cout << "Facet " << b << " corresponds to quadratic boundary "
1703 // << boundary_id[b] << " which contains "
1704 // << this->nboundary_element(boundary_id[b])
1705 // << " element[s] " << std::endl;
1706  }
1707 
1708 
1709  // Loop over bulk elements that are adjacent to quadratic boundary
1710  Vector<double> boundary_zeta(2);
1711  Vector<double> quadratic_coordinates(3);
1712  GeomObject* geom_obj_pt=0;
1713  Vector<double> s_geom_obj(2);
1714  unsigned nel=this->nboundary_element(boundary_id[b]);
1715  for (unsigned e=0;e<nel;e++)
1716  {
1717  // Get pointer to the bulk element that is adjacent to boundary
1718  FiniteElement* bulk_elem_pt=this->boundary_element_pt(boundary_id[b],e);
1719 
1720  // Loop over nodes
1721  unsigned nnod=bulk_elem_pt->nnode();
1722  for (unsigned j=0;j<nnod;j++)
1723  {
1724  Node* nod_pt=bulk_elem_pt->node_pt(j);
1725  if (nod_pt->is_on_boundary(boundary_id[b]))
1726  {
1727  nod_pt->get_coordinates_on_boundary(boundary_id[b],boundary_zeta);
1728 
1729  // Doc it?
1730  if (doc_info.is_doc_enabled())
1731  {
1732  the_file_before << nod_pt->x(0) << " "
1733  << nod_pt->x(1) << " "
1734  << nod_pt->x(2) << " "
1735  << boundary_zeta[0] << " "
1736  << boundary_zeta[1] << " "
1737  << std::endl;
1738  }
1739 
1740  // Find local coordinate in quadratic facet
1741  el_pt->locate_zeta(boundary_zeta,geom_obj_pt,s_geom_obj);
1742  if (geom_obj_pt!=0)
1743  {
1744  geom_obj_pt->position(s_geom_obj,quadratic_coordinates);
1745  nod_pt->x(0)=quadratic_coordinates[0];
1746  nod_pt->x(1)=quadratic_coordinates[1];
1747  nod_pt->x(2)=quadratic_coordinates[2];
1748  }
1749  else
1750  {
1751  // Get ready to bail out below...
1752  success=false;
1753 
1754  std::ostringstream error_message;
1755  error_message
1756  << "Warning: Couldn't find GeomObject during snapping to\n"
1757  << "quadratic surface for boundary " << boundary_id[b]
1758  << ". I'm leaving the node where it was. Will bail out later.\n";
1760  error_message.str(),
1761  "TetgenMesh::snap_to_quadratic_surface()",
1762  OOMPH_EXCEPTION_LOCATION);
1763  if (!the_file_non_snapped.is_open())
1764  {
1765  the_file_non_snapped.open(non_snapped_filename.c_str());
1766  }
1767  the_file_non_snapped << nod_pt->x(0) << " "
1768  << nod_pt->x(1) << " "
1769  << nod_pt->x(2) << " "
1770  << boundary_zeta[0] << " "
1771  << boundary_zeta[1] << " "
1772  << std::endl;
1773  }
1774 
1775  // Doc it?
1776  if (doc_info.is_doc_enabled())
1777  {
1778  the_file_after << nod_pt->x(0) << " "
1779  << nod_pt->x(1) << " "
1780  << nod_pt->x(2) << " "
1781  << boundary_zeta[0] << " "
1782  << boundary_zeta[1] << " "
1783  << std::endl;
1784  }
1785  }
1786  }
1787  }
1788 
1789  // Close doc file
1790  the_file.close();
1791  the_file_before.close();
1792  the_file_after.close();
1793  }
1794 
1795  // Bail out?
1796  if (!success)
1797  {
1798  std::ostringstream error_message;
1799  error_message
1800  << "Warning: Couldn't find GeomObject during snapping to\n"
1801  << "quadratic surface. Bailing out.\n"
1802  << "Nodes that couldn't be snapped are contained in \n"
1803  << "file: " << non_snapped_filename << ".\n"
1804  << "This problem may arise because the switch_normal flag was \n"
1805  << "set wrongly.\n";
1806  throw OomphLibError(
1807  error_message.str(),
1808  OOMPH_CURRENT_FUNCTION,
1809  OOMPH_EXCEPTION_LOCATION);
1810  }
1811 
1812  // Close
1813  if (!the_file_non_snapped.is_open())
1814  {
1815  the_file_non_snapped.close();
1816  }
1817 
1818  // Kill auxiliary FreeStandingFaceElements
1819  for(unsigned e=0;e<nel;e++)
1820  {
1821  delete face_el_pt[e];
1822  face_el_pt[e]=0;
1823  }
1824 
1825  // Kill boundary nodes
1826  unsigned nn=nod_pt.size();
1827  for (unsigned j=0;j<nn;j++)
1828  {
1829  delete nod_pt[j];
1830  }
1831 
1832 }
1833 
1834 
1835 
1836 //========================================================================
1837 /// Non-delaunay split elements that have three faces on a boundary
1838 /// into sons.
1839 //========================================================================
1840 template <class ELEMENT>
1842  TimeStepper* time_stepper_pt)
1843 {
1844 
1845  // Setup boundary lookup scheme if required
1846  if (!Lookup_for_elements_next_boundary_is_setup)
1847  {
1848  setup_boundary_element_info();
1849  }
1850 
1851  // Find out how many nodes we have along each element edge
1852  unsigned nnode_1d=finite_element_pt(0)->nnode_1d();
1853  // Find out the total number of nodes
1854  unsigned nnode = this->finite_element_pt(0)->nnode();
1855 
1856  // At the moment we're only able to deal with nnode_1d=2 or 3.
1857  if ((nnode_1d!=2)&&(nnode_1d!=3))
1858  {
1859  std::ostringstream error_message;
1860  error_message << "Mesh generation from tetgen currently only works\n";
1861  error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
1862  error_message << "for nnode_1d=" << nnode_1d << std::endl;
1863 
1864  throw OomphLibError(error_message.str(),
1865  OOMPH_CURRENT_FUNCTION,
1866  OOMPH_EXCEPTION_LOCATION);
1867  }
1868 
1869  // Map to store how many boundaries elements are located on
1870  std::map<FiniteElement*,unsigned> count;
1871 
1872  // Loop over all boundaries
1873  unsigned nb=this->nboundary();
1874  for (unsigned b=0;b<nb;b++)
1875  {
1876  // Loop over elements next to that boundary
1877  unsigned nel=this->nboundary_element(b);
1878  for (unsigned e=0;e<nel;e++)
1879  {
1880  // Get pointer to element
1881  FiniteElement* el_pt=boundary_element_pt(b,e);
1882 
1883  // Bump up counter
1884  count[el_pt]++;
1885  }
1886  }
1887 
1888  // To avoid having to check the map for all elements (which will
1889  // inflate it to the size of all elements!), move offending elements
1890  // into set
1891  std::set<FiniteElement*> elements_to_be_split;
1892  for (std::map<FiniteElement*,unsigned>::iterator it=count.begin();
1893  it!=count.end();it++)
1894  {
1895  // Get pointer to element: Key
1896  FiniteElement* el_pt=it->first;
1897 
1898  // Does it have to be split?
1899  if (it->second>2)
1900  {
1901  elements_to_be_split.insert(el_pt);
1902  }
1903  }
1904 
1905  // Vector for retained or newly built elements
1906  unsigned nel=this->nelement();
1907  Vector<FiniteElement*> new_el_pt;
1908  new_el_pt.reserve(nel);
1909 
1910  // Now loop over all elements
1911  for (unsigned e=0;e<nel;e++)
1912  {
1913 
1914  // Get pointer to element
1915  FiniteElement* el_pt=this->finite_element_pt(e);
1916 
1917  // Does it have to be split? I.e. is it in the set?
1918  std::set<FiniteElement*>::iterator it=
1919  std::find(elements_to_be_split.begin(),elements_to_be_split.end(),el_pt);
1920 
1921  // It's not in the set, so iterator has reached end
1922  if (it==elements_to_be_split.end())
1923  {
1924  // Carry it across
1925  new_el_pt.push_back(el_pt);
1926  }
1927  // It's in the set of elements to be split
1928  else
1929  {
1930  // Storage for new nodes -- Note: All new nodes are interior and
1931  // therefore cannot be boundary nodes!
1932  Node* node0_pt=0;
1933  Node* node1_pt=0;
1934  Node* node2_pt=0;
1935  Node* node3_pt=0;
1936  Node* node4_pt=0;
1937  Node* node5_pt=0;
1938  Node* node6_pt=0;
1939  Node* node7_pt=0;
1940  Node* node8_pt=0;
1941  Node* node9_pt=0;
1942  Node* node10_pt=0;
1943 
1944  // Create first new element
1945  FiniteElement* el1_pt=new ELEMENT;
1946 
1947  // Copy existing nodes
1948  el1_pt->node_pt(0)=el_pt->node_pt(0);
1949  el1_pt->node_pt(1)=el_pt->node_pt(1);
1950  el1_pt->node_pt(3)=el_pt->node_pt(3);
1951  if (nnode_1d==3)
1952  {
1953  el1_pt->node_pt(4)=el_pt->node_pt(4);
1954  el1_pt->node_pt(6)=el_pt->node_pt(6);
1955  el1_pt->node_pt(9)=el_pt->node_pt(9);
1956  }
1957 
1958  // Create new nodes
1959  // If we have an enriched element then don't need to construct centroid
1960  // node
1961  if(nnode==15)
1962  {
1963  node0_pt = el_pt->node_pt(14);
1964  el1_pt->node_pt(2) = node0_pt;
1965  node5_pt = el1_pt->construct_node(11,time_stepper_pt);
1966  node6_pt = el1_pt->construct_node(13,time_stepper_pt);
1967  node9_pt = el1_pt->construct_node(12,time_stepper_pt);
1968  //Copy others over
1969  el1_pt->node_pt(10) = el_pt->node_pt(10);
1970  }
1971  //If not enriched we do
1972  else
1973  {
1974  node0_pt=el1_pt->construct_node(2,time_stepper_pt);
1975  }
1976  if (nnode_1d==3)
1977  {
1978  node1_pt=el1_pt->construct_boundary_node(5,time_stepper_pt);
1979  node2_pt=el1_pt->construct_boundary_node(7,time_stepper_pt);
1980  node4_pt=el1_pt->construct_boundary_node(8,time_stepper_pt);
1981  }
1982 
1983 
1984  // Create second new element
1985  FiniteElement* el2_pt=new ELEMENT;
1986 
1987  // Copy existing nodes
1988  el2_pt->node_pt(0)=el_pt->node_pt(0);
1989  el2_pt->node_pt(1)=el_pt->node_pt(1);
1990  el2_pt->node_pt(2)=el_pt->node_pt(2);
1991  if (nnode_1d==3)
1992  {
1993  el2_pt->node_pt(4)=el_pt->node_pt(4);
1994  el2_pt->node_pt(5)=el_pt->node_pt(5);
1995  el2_pt->node_pt(7)=el_pt->node_pt(7);
1996  }
1997 
1998  // Create new node
1999  if (nnode_1d==3)
2000  {
2001  node3_pt=el2_pt->construct_boundary_node(8,time_stepper_pt);
2002  }
2003 
2004  // Copy existing new nodes
2005  el2_pt->node_pt(3)=node0_pt;
2006  if (nnode_1d==3)
2007  {
2008  el2_pt->node_pt(6)=node1_pt;
2009  el2_pt->node_pt(9)=node2_pt;
2010  }
2011 
2012  //Copy and constuct other nodes for enriched elements
2013  if(nnode==15)
2014  {
2015  el2_pt->node_pt(10) = node5_pt;
2016  el2_pt->node_pt(11) = el_pt->node_pt(11);
2017  node8_pt = el2_pt->construct_node(12,time_stepper_pt);
2018  node10_pt = el2_pt->construct_node(13,time_stepper_pt);
2019  }
2020 
2021  // Create third new element
2022  FiniteElement* el3_pt=new ELEMENT;
2023 
2024  // Copy existing nodes
2025  el3_pt->node_pt(1)=el_pt->node_pt(1);
2026  el3_pt->node_pt(2)=el_pt->node_pt(2);
2027  el3_pt->node_pt(3)=el_pt->node_pt(3);
2028  if (nnode_1d==3)
2029  {
2030  el3_pt->node_pt(7)=el_pt->node_pt(7);
2031  el3_pt->node_pt(8)=el_pt->node_pt(8);
2032  el3_pt->node_pt(9)=el_pt->node_pt(9);
2033  }
2034 
2035  // Copy existing new nodes
2036  el3_pt->node_pt(0)=node0_pt;
2037  if (nnode_1d==3)
2038  {
2039  el3_pt->node_pt(4)=node2_pt;
2040  el3_pt->node_pt(5)=node3_pt;
2041  el3_pt->node_pt(6)=node4_pt;
2042  }
2043 
2044  //Copy and constuct other nodes for enriched elements
2045  if(nnode==15)
2046  {
2047  el3_pt->node_pt(10) = node6_pt;
2048  el3_pt->node_pt(11) = node10_pt;
2049  node7_pt = el3_pt->construct_node(12,time_stepper_pt);
2050  el3_pt->node_pt(13) = el_pt->node_pt(13);
2051  }
2052 
2053 
2054  // Create fourth new element
2055  FiniteElement* el4_pt=new ELEMENT;
2056 
2057  // Copy existing nodes
2058  el4_pt->node_pt(0)=el_pt->node_pt(0);
2059  el4_pt->node_pt(2)=el_pt->node_pt(2);
2060  el4_pt->node_pt(3)=el_pt->node_pt(3);
2061  if (nnode_1d==3)
2062  {
2063  el4_pt->node_pt(5)=el_pt->node_pt(5);
2064  el4_pt->node_pt(6)=el_pt->node_pt(6);
2065  el4_pt->node_pt(8)=el_pt->node_pt(8);
2066  }
2067 
2068  // Copy existing new nodes
2069  el4_pt->node_pt(1)=node0_pt;
2070  if (nnode_1d==3)
2071  {
2072  el4_pt->node_pt(4)=node1_pt;
2073  el4_pt->node_pt(7)=node3_pt;
2074  el4_pt->node_pt(9)=node4_pt;
2075  }
2076 
2077  //Copy all other nodes
2078  if(nnode==15)
2079  {
2080  el4_pt->node_pt(10) = node9_pt;
2081  el4_pt->node_pt(11) = node8_pt;
2082  el4_pt->node_pt(12) = el_pt->node_pt(12);
2083  el4_pt->node_pt(13) = node7_pt;;
2084  }
2085 
2086 
2087  // Add new elements and nodes
2088  new_el_pt.push_back(el1_pt);
2089  new_el_pt.push_back(el2_pt);
2090  new_el_pt.push_back(el3_pt);
2091  new_el_pt.push_back(el4_pt);
2092 
2093  if(nnode!=15)
2094  {
2095  this->add_node_pt(node0_pt);
2096  }
2097  this->add_node_pt(node1_pt);
2098  this->add_node_pt(node2_pt);
2099  this->add_node_pt(node3_pt);
2100  this->add_node_pt(node4_pt);
2101  if(nnode==15)
2102  {
2103  this->add_node_pt(node5_pt);
2104  this->add_node_pt(node6_pt);
2105  this->add_node_pt(node7_pt);
2106  this->add_node_pt(node8_pt);
2107  this->add_node_pt(node9_pt);
2108  }
2109 
2110  // Set nodal positions
2111  for (unsigned i=0;i<3;i++)
2112  {
2113  //Only bother to set centroid if does not already exist
2114  if(nnode!=15)
2115  {
2116  node0_pt->x(i)=0.25*(el_pt->node_pt(0)->x(i)+
2117  el_pt->node_pt(1)->x(i)+
2118  el_pt->node_pt(2)->x(i)+
2119  el_pt->node_pt(3)->x(i));
2120  }
2121 
2122  if (nnode_1d==3)
2123  {
2124  node1_pt->x(i)=0.5*(el_pt->node_pt(0)->x(i)+node0_pt->x(i));
2125  node2_pt->x(i)=0.5*(el_pt->node_pt(1)->x(i)+node0_pt->x(i));
2126  node3_pt->x(i)=0.5*(el_pt->node_pt(2)->x(i)+node0_pt->x(i));
2127  node4_pt->x(i)=0.5*(el_pt->node_pt(3)->x(i)+node0_pt->x(i));
2128  }
2129  }
2130 
2131 
2132  //Construct the four interior nodes if needed
2133  //and add to the list
2134  if(nnode==15)
2135  {
2136  //Set the positions of the newly created mid-face nodes
2137  //New Node 5 lies in the plane between original nodes 0 1 centroid
2138  for(unsigned i=0;i<3;++i)
2139  {
2140  node5_pt->x(i) =
2141  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i)
2142  + el_pt->node_pt(14)->x(i))/3.0;
2143  }
2144 
2145  //New Node 6 lies in the plane between original nodes 1 3 centroid
2146  for(unsigned i=0;i<3;++i)
2147  {
2148  node6_pt->x(i) =
2149  (el_pt->node_pt(1)->x(i) + el_pt->node_pt(3)->x(i)
2150  + el_pt->node_pt(14)->x(i))/3.0;
2151  }
2152 
2153  //New Node 7 lies in the plane between original nodes 2 3 centroid
2154  for(unsigned i=0;i<3;++i)
2155  {
2156  node7_pt->x(i) =
2157  (el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i)
2158  + el_pt->node_pt(14)->x(i))/3.0;
2159  }
2160 
2161  //New Node 8 lies in the plane between original nodes 0 2 centroid
2162  for(unsigned i=0;i<3;++i)
2163  {
2164  node8_pt->x(i) =
2165  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(2)->x(i)
2166  + el_pt->node_pt(14)->x(i))/3.0;
2167  }
2168 
2169  //New Node 9 lies in the plane between original nodes 0 3 centroid
2170  for(unsigned i=0;i<3;++i)
2171  {
2172  node9_pt->x(i) =
2173  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(3)->x(i)
2174  + el_pt->node_pt(14)->x(i))/3.0;
2175  }
2176 
2177  //New Node 10 lies in the plane between original nodes 1 2 centroid
2178  for(unsigned i=0;i<3;++i)
2179  {
2180  node10_pt->x(i) =
2181  (el_pt->node_pt(1)->x(i) + el_pt->node_pt(2)->x(i)
2182  + el_pt->node_pt(14)->x(i))/3.0;
2183  }
2184 
2185  //Now create the new centroid nodes
2186 
2187  //First element
2188  Node* temp_nod_pt = el1_pt->construct_node(14,time_stepper_pt);
2189  for(unsigned i=0;i<3;++i)
2190  {
2191  double av_pos = 0.0;
2192  for(unsigned j=0;j<4;j++)
2193  {
2194  av_pos += el1_pt->node_pt(j)->x(i);
2195  }
2196 
2197  temp_nod_pt->x(i) = 0.25*av_pos;
2198  }
2199 
2200  this->add_node_pt(temp_nod_pt);
2201 
2202  //Second element
2203  temp_nod_pt = el2_pt->construct_node(14,time_stepper_pt);
2204  for(unsigned i=0;i<3;++i)
2205  {
2206  double av_pos = 0.0;
2207  for(unsigned j=0;j<4;j++)
2208  {
2209  av_pos += el2_pt->node_pt(j)->x(i);
2210  }
2211  temp_nod_pt->x(i) = 0.25*av_pos;
2212  }
2213  this->add_node_pt(temp_nod_pt);
2214 
2215  //Third element
2216  temp_nod_pt = el3_pt->construct_node(14,time_stepper_pt);
2217  for(unsigned i=0;i<3;++i)
2218  {
2219  double av_pos = 0.0;
2220  for(unsigned j=0;j<4;j++)
2221  {
2222  av_pos += el3_pt->node_pt(j)->x(i);
2223  }
2224  temp_nod_pt->x(i) = 0.25*av_pos;
2225  }
2226  this->add_node_pt(temp_nod_pt);
2227 
2228  //Fourth element
2229  temp_nod_pt = el4_pt->construct_node(14,time_stepper_pt);
2230  for(unsigned i=0;i<3;++i)
2231  {
2232  double av_pos = 0.0;
2233  for(unsigned j=0;j<4;j++)
2234  {
2235  av_pos += el4_pt->node_pt(j)->x(i);
2236  }
2237  temp_nod_pt->x(i) = 0.25*av_pos;
2238  }
2239  this->add_node_pt(temp_nod_pt);
2240  }
2241 
2242 // std::ofstream junk("junk.dat");
2243 // el_pt->output(junk);
2244 // el1_pt->output(junk);
2245 // el2_pt->output(junk);
2246 // el3_pt->output(junk);
2247 // el4_pt->output(junk);
2248 // junk.close();
2249 
2250  // Kill old element
2251  delete el_pt;
2252 
2253  }
2254  }
2255 
2256  // Flush element storage
2257  Element_pt.clear();
2258 
2259  // Copy across
2260  nel=new_el_pt.size();
2261  Element_pt.resize(nel);
2262  for (unsigned e=0;e<nel;e++)
2263  {
2264  Element_pt[e]=new_el_pt[e];
2265  }
2266 
2267  // Setup boundary lookup scheme again
2268  setup_boundary_element_info();
2269 
2270 }
2271 
2272 }
2273 #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
void split_elements_in_corners(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Non-Delaunay split elements that have three faces on a boundary into sons. Timestepper species timest...
void build_from_scaffold(TimeStepper *time_stepper_pt, const bool &use_attributes)
Build mesh from scaffold.
bool is_doc_enabled() const
Are we documenting?
Information for documentation of results: Directory and file number to enable output in the form RESL...
cstr elem_len * i
Definition: cfortran.h:607
void snap_to_quadratic_surface(const Vector< unsigned > &boundary_id, const std::string &quadratic_surface_file_name, const bool &switch_normal, DocInfo &doc_info)
Snap boundaries specified by the IDs listed in boundary_id to a quadratric surface, specified in the file quadratic_surface_file_name. This is usually used with vmtk-based meshes for which oomph-lib's xda to poly conversion code produces the files "quadratic_fsi_boundary.dat" and "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI boundary (for the fluid and the solid) and the quadratic representation of the outer boundary of the solid. When used with these files, the flag switch_normal should be set to true when calling the function for the outer boundary of the solid. The DocInfo object can be used to label optional output files. (Uses directory and label).
std::string & label()
String used (e.g.) for labeling output files.
A general Finite Element class.
Definition: elements.h:1271
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary, so return false. This will be overloaded by BoundaryNodes.
Definition: nodes.h:1293
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition: elements.h:2408
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
e
Definition: cfortran.h:575
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Return the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interface provides run-time error checking.
Definition: nodes.cc:2287
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
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object. ...
Definition: elements.h:2380
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh...
Definition: nodes.h:67
static char t char * s
Definition: cfortran.h:572
void deep_copy_of_tetgenio(tetgenio *const &input_pt, tetgenio *&output_pt)
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
std::string directory() const
Output directory.
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual) ...
Definition: elements.h:1795
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual void setup_boundary_coordinates_generic(const unsigned &b, const bool &switch_normal, std::ofstream &outfile)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219