tet_mesh.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Common base class for all Tet Meshes
27 #ifndef OOMPH_TET_MESH_HEADER
28 #define OOMPH_TET_MESH_HEADER
29 
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32 #include <oomph-lib-config.h>
33 #endif
34 
35 // oomph-lib includes
36 #include "Vector.h"
37 #include "nodes.h"
38 #include "matrices.h"
39 #include "mesh.h"
40 #include "geom_obj_with_boundary.h"
41 
42 namespace oomph
43 {
44  //=======================================================================
45  /// Vertex for Tet mesh generation. Can lie on multiple boundaries
46  /// (identified via one-based enumeration!) and can have intrinisic
47  /// coordinates in a DiskLikeGeomObjectWithBoundaries.
48  //=======================================================================
50  {
51  public:
52  /// Only friends can set boundary ID -- the facet is my only friend!
53  friend class TetMeshFacet;
54 
55  /// Constructor: Pass coordinates (length 3!)
57  {
58 #ifdef PARANOID
59  if (X.size() != 3)
60  {
61  std::ostringstream error_stream;
62  error_stream << "TetMeshVertex should only be used in 3D!\n"
63  << "Your Vector of coordinates, contains data for "
64  << x.size() << "-dimensional coordinates." << std::endl;
65  throw OomphLibError(
66  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
67  }
68 #endif
69  }
70 
71  // Constructor: Take coordinates from a node
72  TetMeshVertex(Node* const& node_pt)
73  {
74  const unsigned n_dim = node_pt->ndim();
75 #ifdef PARANOID
76  if (n_dim != 3)
77  {
78  std::ostringstream error_stream;
79  error_stream << "TetMeshVertex should only be used in 3D!\n"
80  << "Your Node contains data for " << n_dim
81  << "-dimensional coordinates." << std::endl;
82  throw OomphLibError(
83  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
84  }
85 #endif
86 
87  // Copy over the positions from the node
88  X.resize(n_dim);
89  for (unsigned i = 0; i < n_dim; ++i)
90  {
91  X[i] = node_pt->x(i);
92  }
93  }
94 
95 
96  /// Set intrinisic coordinates in GeomObject
98  {
99 #ifdef PARANOID
100  if (zeta.size() != 2)
101  {
102  std::ostringstream error_stream;
103  error_stream
104  << "TetMeshVertex should only be used in 3D!\n"
105  << "Your Vector of intrinisic coordinates, contains data for "
106  << zeta.size() << "-dimensional coordinates but should be 2!"
107  << std::endl;
108  throw OomphLibError(
109  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
110  }
111 #endif
112 
113  Zeta_in_geom_object = zeta;
114  }
115 
116  /// Get intrinisic coordinates in GeomObject (returns zero sized
117  /// vector if no such coordinates have been specified)
119  {
120  return Zeta_in_geom_object;
121  }
122 
123  /// i-th coordinate
124  double x(const unsigned& i) const
125  {
126  return X[i];
127  }
128 
129  /// First (of possibly multiple) one-based boundary id
130  unsigned one_based_boundary_id() const
131  {
132  if (One_based_boundary_id.size() == 0)
133  {
134  return 0;
135  }
136  return *(One_based_boundary_id.begin());
137  }
138 
139  private:
140  /// Set of (one-based!) boundary IDs this vertex lives on
141  void set_one_based_boundary_id(const unsigned& id)
142  {
143  One_based_boundary_id.insert(id);
144  }
145 
146  /// Coordinate vector
148 
149  /// Set of (one-based!) boundary IDs this vertex lives on
150  std::set<unsigned> One_based_boundary_id;
151 
152  /// Intrinisic coordinates in GeomObject (zero sized if there
153  /// isn't one.
155  };
156 
157  /// /////////////////////////////////////////////////////////////////////
158  /// /////////////////////////////////////////////////////////////////////
159  /// /////////////////////////////////////////////////////////////////////
160 
161 
162  //=======================================================================
163  /// Facet for Tet mesh generation. Can lie on boundary
164  /// (identified via one-based enumeration!) and can have
165  /// GeomObject associated with those boundaries.
166  //=======================================================================
168  {
169  public:
170  /// Constructor: Specify number of vertices
171  TetMeshFacet(const unsigned& nvertex)
172  : One_based_boundary_id(0), // Initialision implies not on any boundary
174  0) // Initialisation implies
175  // not embedded in any region
176  {
177  Vertex_pt.resize(nvertex, 0);
178  }
179 
180  /// Number of vertices
181  unsigned nvertex() const
182  {
183  return Vertex_pt.size();
184  }
185 
186  /// Pointer to j-th vertex (const)
187  TetMeshVertex* vertex_pt(const unsigned& j) const
188  {
189  return Vertex_pt[j];
190  }
191 
192  /// Set pointer to j-th vertex and pass one-based boundary id that
193  /// may already have been set for this facet.
194  void set_vertex_pt(const unsigned& j, TetMeshVertex* vertex_pt)
195  {
196  Vertex_pt[j] = vertex_pt;
197  // If not set yet, this is harmless since it simply over-writes
198  // the dummy value in vertex
199  TetMeshVertex* v_pt = Vertex_pt[j];
200  if (v_pt != 0)
201  {
203  }
204  }
205 
206  /// Constant access to (one-based!) boundary IDs this facet lives on
207  unsigned one_based_boundary_id() const
208  {
209  return One_based_boundary_id;
210  }
211 
212  /// Set (one-based!) boundary IDs this facet lives on. Passed to any
213  /// existing vertices and to any future ones
214  void set_one_based_boundary_id(const unsigned& one_based_id)
215  {
216  One_based_boundary_id = one_based_id;
217  unsigned nv = Vertex_pt.size();
218  for (unsigned j = 0; j < nv; j++)
219  {
220  TetMeshVertex* v_pt = Vertex_pt[j];
221  if (v_pt != 0)
222  {
223  v_pt->set_one_based_boundary_id(one_based_id);
224  }
225  }
226  }
227 
228  /// Set (one-based!) region ID this facet is adjacent to.
229  /// Specification of zero argument is harmless as it indicates that
230  /// that facet is not adjacent to any region.
231  void set_one_based_adjacent_region_id(const unsigned& one_based_id)
232  {
233  One_based_adjacent_region_id.insert(one_based_id);
234  }
235 
236  /// Return set of (one-based!) region IDs this facet is adjacent to
237  std::set<unsigned> one_based_adjacent_region_id() const
238  {
240  }
241 
242  /// Boolean indicating that facet is embedded in a specified region
244  {
246  }
247 
248  /// Facet is to be embedded in specified one-based region
250  const unsigned& one_based_region_id)
251  {
252  One_based_region_id_that_facet_is_embedded_in = one_based_region_id;
253  }
254 
255  /// Which (one-based) region is facet embedded in (if zero, it's not
256  /// embedded in any region)
258  {
260  }
261 
262  /// Output
263  void output(std::ostream& outfile) const
264  {
265  unsigned n = Vertex_pt.size();
266  outfile << "ZONE I=" << n + 1 << std::endl;
267  for (unsigned j = 0; j < n; j++)
268  {
269  outfile << Vertex_pt[j]->x(0) << " " << Vertex_pt[j]->x(1) << " "
270  << Vertex_pt[j]->x(2) << " " << One_based_boundary_id
271  << std::endl;
272  }
273  outfile << Vertex_pt[0]->x(0) << " " << Vertex_pt[0]->x(1) << " "
274  << Vertex_pt[0]->x(2) << " " << One_based_boundary_id
275  << std::endl;
276  }
277 
278  private:
279  /// Pointer to vertices
281 
282  /// (One-based!) boundary IDs this facet lives on
284 
285  /// Set of one-based adjacent region ids; no adjacent region if
286  /// it's zero.
287  std::set<unsigned> One_based_adjacent_region_id;
288 
289 
290  /// Facet is to be embedded in specified one-based region.
291  /// Defaults to zero, indicating that its not embedded.
293  };
294 
295 
296  /// /////////////////////////////////////////////////////////////////////
297  /// /////////////////////////////////////////////////////////////////////
298  /// /////////////////////////////////////////////////////////////////////
299 
300 
301  //========================================================================
302  /// Base class for tet mesh boundary defined by polygonal
303  /// planar facets
304  //========================================================================
306  {
307  public:
308  /// Constructor:
312  {
313  }
314 
315  /// Empty destructor
317 
318  /// Number of vertices
319  unsigned nvertex() const
320  {
321  return Vertex_pt.size();
322  }
323 
324  /// Number of facets
325  unsigned nfacet() const
326  {
327  return Facet_pt.size();
328  }
329 
330  /// One-based boundary id of j-th facet
331  unsigned one_based_facet_boundary_id(const unsigned& j) const
332  {
333  return Facet_pt[j]->one_based_boundary_id();
334  }
335 
336  /// First (of possibly multiple) one-based boundary id of j-th vertex
337  unsigned one_based_vertex_boundary_id(const unsigned& j) const
338  {
339  return Vertex_pt[j]->one_based_boundary_id();
340  }
341 
342  /// i-th coordinate of j-th vertex
343  double vertex_coordinate(const unsigned& j, const unsigned& i) const
344  {
345  return Vertex_pt[j]->x(i);
346  }
347 
348  /// Number of vertices defining the j-th facet
349  unsigned nvertex_on_facet(const unsigned& j) const
350  {
351  return Facet_pt[j]->nvertex();
352  }
353 
354  /// Test whether boundary can be split in tetgen
356  {
358  }
359 
360  /// Test whether boundaries can be split in tetgen
362  {
364  }
365 
366  /// Test whether boundaries can be split in tetgen
368  {
370  }
371 
372  /// Pointer to j-th facet
373  TetMeshFacet* facet_pt(const unsigned& j) const
374  {
375  return Facet_pt[j];
376  }
377 
378  /// Pointer to j-th vertex
379  TetMeshVertex* vertex_pt(const unsigned& j) const
380  {
381  return Vertex_pt[j];
382  }
383 
384  /// Access to GeomObject with boundaries associated with this
385  /// surface (Null if there isn't one!)
387  {
389  }
390 
391  /// Output
392  void output(std::ostream& outfile) const
393  {
394  unsigned n = Facet_pt.size();
395  for (unsigned j = 0; j < n; j++)
396  {
397  Facet_pt[j]->output(outfile);
398  }
399  }
400 
401 
402  /// Output
403  void output(const std::string& filename) const
404  {
405  std::ofstream outfile;
406  outfile.open(filename.c_str());
407  output(outfile);
408  outfile.close();
409  }
410 
411 
412  //========================================================
413  /// Outputs the faceted surface into a specified file
414  /// in the Paraview format for viewing in Paraview.
415  /// Make sure to output the file with a .vtu extension.
416  /// (Not particularly optimised)
417  //========================================================
418  void output_paraview(std::ostream& outfile) const
419  {
420  // Create storage for vertices, connectivity, offsets and types. These are
421  // all outputted into the .vtu file.
422  std::set<TetMeshVertex*> vertices;
423  Vector<unsigned> connectivity;
424  Vector<unsigned> offsets;
425  Vector<unsigned> types;
426 
427  // Storage for the number of vertices on a facet
428  unsigned n_vertices = 0;
429 
430  // Add every vertex to a set
431  unsigned n_facets = Facet_pt.size();
432  for (unsigned i = 0; i < n_facets; i++)
433  {
434  // Loop over all the vertices of the ith facet and add them to the set
435  // of vertices.
436  n_vertices = Facet_pt[i]->nvertex();
437  for (unsigned j = 0; j < n_vertices; j++)
438  {
439  vertices.insert(Facet_pt[i]->vertex_pt(j));
440  }
441  }
442 
443  // Store the offset of the current facet
444  unsigned current_offset = 0;
445 
446  // This iterator is used for finding the index of a chosen vertex in
447  // 'vertices'
448  std::set<TetMeshVertex*>::iterator iterator;
449 
450  // Get the offset, types and connectivity
451  for (const auto& individual_facet_pt : Facet_pt)
452  {
453  // Get the type of the current facet
454  n_vertices = individual_facet_pt->nvertex();
455 
456  // The types of cell (polygon, tetra, etc) in the VTK/U file format are
457  // enumerated. The type of each facet is inputted using its
458  // corresponding enumeration into the 'types' section of the .vtu file.
459 
460  // The relevant types of cells and enumerations are listed below.
461  // VTK_TRIANGLE = 5
462  // VTK_QUAD = 9
463  // VTK_POLYGON = 7
464 
465  // The full list of types with enumerations are listed on the following
466  // site: https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
467 
468  // Append the value corresponding to the current facet's type.
469  if (n_vertices == 3)
470  {
471  types.push_back(5);
472  }
473  else if (n_vertices == 4)
474  {
475  types.push_back(9);
476  }
477  else
478  {
479  types.push_back(7);
480  }
481 
482  // Find the connectivity of this facet, the 'connectivity' vector
483  // consists of indices specifying the vertices that connect to create
484  // facets.
485  for (unsigned j = 0; j < n_vertices; j++)
486  {
487  iterator = vertices.find(individual_facet_pt->vertex_pt(j));
488  connectivity.push_back(std::distance(vertices.begin(), iterator));
489  }
490 
491  // The vector 'offset' specifies which consecutive elements of
492  // 'connectivity' refer to a single facet.
493  current_offset += n_vertices;
494  offsets.push_back(current_offset);
495  }
496 
497  //========================================================================
498  // Output to file
499  //========================================================================
500 
501  // File Declaration
502  //-----------------
503 
504  // Insert the necessary lines plus the number of vertices and cells/facets
505  // The number of facets is equal to the size of 'offsets'
506  outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
507  "byte_order=\"BigEndian\">\n"
508  << "<UnstructuredGrid>\n"
509  << "<Piece NumberOfPoints=\"" << vertices.size()
510  << "\" NumberOfCells=\"" << offsets.size() << "\">" << std::endl;
511 
512  // Vertices
513  //---------
514 
515  outfile << "<Points>\n"
516  << "<DataArray type=\"Float32\" NumberOfComponents=\"3\" "
517  "format=\"ascii\">"
518  << std::endl;
519 
520  // Output the coordinates of every distinct vertex
521  for (const auto& vertex : vertices)
522  {
523  outfile << vertex->x(0) << " " << vertex->x(1) << " " << vertex->x(2)
524  << std::endl;
525  }
526 
527  outfile << "</DataArray>\n</Points>" << std::endl;
528 
529  // Connectivity
530  //-------------
531 
532  outfile << "<Cells>\n<DataArray type=\"Int32\" "
533  "Name=\"connectivity\" format=\"ascii\">"
534  << std::endl;
535 
536  // Output the connectivity
537  for (const auto& vertex_index : connectivity)
538  {
539  outfile << vertex_index << " ";
540  }
541  outfile << std::endl;
542 
543  outfile << "</DataArray>" << std::endl;
544 
545  // Offsets
546  //--------
547 
548  outfile << "<DataArray type=\"Int32\" Name=\"offsets\" "
549  "format=\"ascii\">"
550  << std::endl;
551 
552  // Output the offsets
553  for (const auto& offset : offsets)
554  {
555  outfile << offset << " ";
556  }
557  outfile << std::endl;
558 
559  outfile << "</DataArray>" << std::endl;
560 
561  // Types
562  //------
563 
564  outfile << "<DataArray type=\"Int32\" Name=\"types\" "
565  "format=\"ascii\">"
566  << std::endl;
567 
568  // Output the types
569  for (const auto& type : types)
570  {
571  outfile << type << " ";
572  }
573  outfile << std::endl;
574 
575  outfile << "</DataArray>\n</Cells>" << std::endl;
576 
577  // File closure
578  //-------------
579 
580  outfile << "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
581  }
582 
583 
584  //========================================================
585  /// Outputs the faceted surface into a file with the
586  /// specified name in the Paraview format.
587  /// (Not particularly optimised)
588  //========================================================
589  void output_paraview(const std::string& filename) const
590  {
591  std::ofstream outfile;
592  outfile.open(filename.c_str());
593  output_paraview(outfile);
594  outfile.close();
595  }
596 
597  /// Virtual function that specifies the variation of the
598  /// zeta coordinates in the GeomObject along the
599  /// boundary connecting vertices 0 and 1 in
600  /// facet facet_id. Default implementation: Linear interpolation
601  /// between edges. zeta_boundary=0.0: we're on vertex 0;
602  /// zeta_boundary=1.0: we're on vertex 1.
603  virtual void boundary_zeta01(const unsigned& facet_id,
604  const double& zeta_boundary,
605  Vector<double>& zeta)
606  {
607  Vector<Vector<double>> zeta_vertex(2);
608  zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(0)->zeta_in_geom_object();
609  zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(1)->zeta_in_geom_object();
610  zeta[0] = zeta_vertex[0][0] +
611  (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
612  zeta[1] = zeta_vertex[0][1] +
613  (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
614  }
615 
616  /// Virtual function that specifies the variation of the
617  /// zeta coordinates in the GeomObject along the
618  /// boundary connecting vertices 1 and 2 in
619  /// facet facet_id. Default implementation: Linear interpolation
620  /// between edges. zeta_boundary=0.0: we're on vertex 1;
621  /// zeta_boundary=1.0: we're on vertex 2.
622  virtual void boundary_zeta12(const unsigned& facet_id,
623  const double& zeta_boundary,
624  Vector<double>& zeta)
625  {
626  Vector<Vector<double>> zeta_vertex(2);
627  zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(1)->zeta_in_geom_object();
628  zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(2)->zeta_in_geom_object();
629  zeta[0] = zeta_vertex[0][0] +
630  (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
631  zeta[1] = zeta_vertex[0][1] +
632  (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
633  }
634 
635  /// Virtual function that specifies the variation of the
636  /// zeta coordinates in the GeomObject along the
637  /// boundary connecting vertices 2 and 0 in
638  /// facet facet_id. Default implementation: Linear interpolation
639  /// between edges. zeta_boundary=0.0: we're on vertex 2;
640  /// zeta_boundary=1.0: we're on vertex 0.
641  virtual void boundary_zeta20(const unsigned& facet_id,
642  const double& zeta_boundary,
643  Vector<double>& zeta)
644  {
645  Vector<Vector<double>> zeta_vertex(2);
646  zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(2)->zeta_in_geom_object();
647  zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(0)->zeta_in_geom_object();
648  zeta[0] = zeta_vertex[0][0] +
649  (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
650  zeta[1] = zeta_vertex[0][1] +
651  (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
652  }
653 
654 
655  /// Facet connectivity: vertex_index[j] is the index of the
656  /// j-th vertex (in the Vertex_pt vector) in facet f. Bit of an obscure
657  /// functionality that's only needed for setup tetgen_io.
659  {
660  if (Facet_vertex_index_in_tetgen.size() != nfacet())
661  {
663  }
665  }
666 
667  protected:
668  /// Vector pointers to vertices
670 
671  /// Vector of pointers to facets
673 
674  /// Boolean to indicate whether extra vertices can be added
675  /// on the boundary in tetgen
677 
678  /// Facet connectivity: Facet_vertex_index[f][j] is the index of the
679  /// j-th vertex (in the Vertex_pt vector) in facet f.
681 
682  /// GeomObject with boundaries associated with this surface
684 
685 
686  private:
687  /// Setup facet connectivity for tetgen
689  {
690  unsigned nv_overall = Vertex_pt.size();
691  unsigned nf = nfacet();
692  Facet_vertex_index_in_tetgen.resize(nf);
693  for (unsigned f = 0; f < nf; f++)
694  {
695  unsigned nv = Facet_pt[f]->nvertex();
696  Facet_vertex_index_in_tetgen[f].resize(nv);
697  for (unsigned v = 0; v < nv; v++)
698  {
699  TetMeshVertex* my_vertex_pt = Facet_pt[f]->vertex_pt(v);
700  for (unsigned j = 0; j < nv_overall; j++)
701  {
702  if (my_vertex_pt == Vertex_pt[j])
703  {
705  break;
706  }
707  }
708  }
709  }
710  }
711  };
712 
713 
714  /// /////////////////////////////////////////////////////////////////////
715  /// /////////////////////////////////////////////////////////////////////
716  /// /////////////////////////////////////////////////////////////////////
717 
718 
719  //========================================================================
720  /// Base class for closed tet mesh boundary bounded by polygonal
721  /// planar facets
722  //========================================================================
724  {
725  public:
726  /// Constructor:
729  {
730  }
731 
732  /// Empty destructor
734 
735  /// Declare closed surface to represent hole for gmsh
737  {
739  }
740 
741  /// Declare closed surface NOT to represent hole for gmsh
743  {
745  }
746 
747  /// Does closed surface represent hole for gmsh?
749  {
751  }
752 
753  /// i=th coordinate of the j-th internal point for tetgen
754  const double& internal_point_for_tetgen(const unsigned& j,
755  const unsigned& i) const
756  {
757  return (Internal_point_for_tetgen[j].first)[i];
758  }
759 
760  /// Specify coordinate of hole for tetgen
761  void set_hole_for_tetgen(const Vector<double>& hole_point)
762  {
763  Internal_point_for_tetgen.push_back(std::make_pair(hole_point, -1));
764  }
765 
766  /// Specify a region; pass (zero-based) region ID and coordinate
767  /// of point in region for tetgen
768  void set_region_for_tetgen(const unsigned& region_id,
769  const Vector<double>& region_point)
770  {
771  Internal_point_for_tetgen.push_back(
772  std::make_pair(region_point, region_id));
773  }
774 
775  /// Number of internal points (identifying either regions or holes)
776  /// for tetgen
778  {
779  return Internal_point_for_tetgen.size();
780  }
781 
782  /// Return the (zero-based) region ID of j-th internal point for
783  /// tetgen. Negative if it's actually a hole.
784  const int& region_id_for_tetgen(const unsigned& j) const
785  {
786  return Internal_point_for_tetgen[j].second;
787  }
788 
789 
790  /// Is j-th internal point for tetgen associated with a hole?
792  {
793  return (Internal_point_for_tetgen[j].second < 0);
794  }
795 
796  /// Is j-th internal point for tetgen associated with a region?
798  {
799  return (Internal_point_for_tetgen[j].second >= 0);
800  }
801 
802 
803  private:
804  /// Storage for internal points for tetgen. Stores pair of:
805  /// -- Vector containing coordinates of internal point
806  /// -- region ID (negative if it's a hole)
808 
809  /// Does closed surface represent hole for gmsh?
811  };
812 
813  /// ////////////////////////////////////////////////////////////////////
814  /// ////////////////////////////////////////////////////////////////////
815  /// ////////////////////////////////////////////////////////////////////
816 
817 
820  {
821  public:
822  // Constructor, which requires node, connectivity and boundary information
824  Vector<Node*> const& vertex_node_pt,
825  Vector<Vector<unsigned>> const& facet_connectivity,
826  Vector<unsigned> const& facet_boundary_id);
827 
828  // Destructor
830  };
831 
832 
833  /// ////////////////////////////////////////////////////////////////////
834  /// ////////////////////////////////////////////////////////////////////
835  /// ////////////////////////////////////////////////////////////////////
836 
837 
838  /// ////////////////////////////////////////////////////////////////////
839  /// ////////////////////////////////////////////////////////////////////
840  /// ////////////////////////////////////////////////////////////////////
841 
842 
843  //================================================================
844  /// Base class for tet meshes (meshes made of 3D tet elements).
845  //================================================================
846  class TetMeshBase : public virtual Mesh
847  {
848  public:
849  /// Constructor
851 
852  /// Broken copy constructor
853  TetMeshBase(const TetMeshBase& node) = delete;
854 
855  /// Broken assignment operator
856  void operator=(const TetMeshBase&) = delete;
857 
858  /// Destructor (empty)
859  virtual ~TetMeshBase() {}
860 
861  /// Global static data that specifies the permitted
862  /// error in the setup of the boundary coordinates
864 
865  /// Assess mesh quality: Ratio of max. edge length to min. height,
866  /// so if it's very large it's BAAAAAD.
867  void assess_mesh_quality(std::ofstream& some_file);
868 
869  /// Setup boundary coordinate on boundary b which is
870  /// assumed to be planar. Boundary coordinates are the
871  /// x-y coordinates in the plane of that boundary, with the
872  /// x-axis along the line from the (lexicographically)
873  /// "lower left" to the "upper right" node. The y axis
874  /// is obtained by taking the cross-product of the positive
875  /// x direction with the outer unit normal computed by
876  /// the face elements (or its negative if switch_normal is set
877  /// to true). Doc faces in output file (if it's open).
878  ///
879  /// Note 1: Setup of boundary coordinates is not done if the boundary in
880  /// question turns out to be nonplanar.
881  ///
882  /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
883  /// also
884  /// store the boundary coordinates of its vertices. They are needed
885  /// to interpolated intrinsic coordinates of an associated
886  /// GeomObject (if any) into the interior.
887  template<class ELEMENT>
888  void setup_boundary_coordinates(const unsigned& b)
889  {
890  // Dummy file
891  std::ofstream some_file;
892 
893  // Don't switch the normal
894  bool switch_normal = false;
895  this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, some_file);
896  }
897 
898  /// Setup boundary coordinate on boundary b which is
899  /// assumed to be planar. Boundary coordinates are the
900  /// x-y coordinates in the plane of that boundary, with the
901  /// x-axis along the line from the (lexicographically)
902  /// "lower left" to the "upper right" node. The y axis
903  /// is obtained by taking the cross-product of the positive
904  /// x direction with the outer unit normal computed by
905  /// the face elements (or its negative if switch_normal is set
906  /// to true). Doc faces in output file (if it's open).
907  ///
908  /// Note 1: Setup of boundary coordinates is not done if the boundary in
909  /// question turns out to be nonplanar.
910  ///
911  /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
912  /// also
913  /// store the boundary coordinates of its vertices. They are needed
914  /// to interpolated intrinsic coordinates of an associated
915  /// GeomObject (if any) into the interior.
916  /// Final boolean argument allows switching of the direction of the outer
917  /// unit normal.
918  template<class ELEMENT>
919  void setup_boundary_coordinates(const unsigned& b,
920  const bool& switch_normal)
921  {
922  // Dummy file
923  std::ofstream some_file;
924 
925  this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, some_file);
926  }
927 
928 
929  /// Setup boundary coordinate on boundary b which is
930  /// assumed to be planar. Boundary coordinates are the
931  /// x-y coordinates in the plane of that boundary, with the
932  /// x-axis along the line from the (lexicographically)
933  /// "lower left" to the "upper right" node. The y axis
934  /// is obtained by taking the cross-product of the positive
935  /// x direction with the outer unit normal computed by
936  /// the face elements (or its negative if switch_normal is set
937  /// to true). Doc faces in output file (if it's open).
938  ///
939  /// Note 1: Setup of boundary coordinates is not done if the boundary in
940  /// question turns out to be nonplanar.
941  ///
942  /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
943  /// also
944  /// store the boundary coordinates of its vertices. They are needed
945  /// to interpolated intrinsic coordinates of an associated
946  /// GeomObject (if any) into the interior.
947  /// Boolean argument allows switching of the direction of the outer
948  /// unit normal. Output file for doc.
949  template<class ELEMENT>
950  void setup_boundary_coordinates(const unsigned& b,
951  const bool& switch_normal,
952  std::ofstream& outfile);
953 
954  /// Setup boundary coordinate on boundary b which is
955  /// assumed to be planar. Boundary coordinates are the
956  /// x-y coordinates in the plane of that boundary, with the
957  /// x-axis along the line from the (lexicographically)
958  /// "lower left" to the "upper right" node. The y axis
959  /// is obtained by taking the cross-product of the positive
960  /// x direction with the outer unit normal computed by
961  /// the face elements (or its negative if switch_normal is set
962  /// to true). Doc faces in output file (if it's open).
963  ///
964  /// Note 1: Setup of boundary coordinates is not done if the boundary in
965  /// question turns out to be nonplanar.
966  ///
967  /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
968  /// also
969  /// store the boundary coordinates of its vertices. They are needed
970  /// to interpolated intrinsic coordinates of an associated
971  /// GeomObject (if any) into the interior.
972  /// Output file for doc.
973  template<class ELEMENT>
974  void setup_boundary_coordinates(const unsigned& b, std::ofstream& outfile)
975  {
976  // Don't switch the normal
977  bool switch_normal = false;
978  this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, outfile);
979  }
980 
981 
982  /// Return the number of elements adjacent to boundary b in region r
983  inline unsigned nboundary_element_in_region(const unsigned& b,
984  const unsigned& r) const
985  {
986  // Need to use a constant iterator here to keep the function "const"
987  // Return an iterator to the appropriate entry, if we find it
988  std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
989  Boundary_region_element_pt[b].find(r);
990  if (it != Boundary_region_element_pt[b].end())
991  {
992  return (it->second).size();
993  }
994  // Otherwise there are no elements adjacent to boundary b in the region r
995  else
996  {
997  return 0;
998  }
999  }
1000 
1001  /// Return pointer to the e-th element adjacent to boundary b in region r
1003  const unsigned& r,
1004  const unsigned& e) const
1005  {
1006  // Use a constant iterator here to keep function "const" overall
1007  std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
1008  Boundary_region_element_pt[b].find(r);
1009  if (it != Boundary_region_element_pt[b].end())
1010  {
1011  return (it->second)[e];
1012  }
1013  else
1014  {
1015  return 0;
1016  }
1017  }
1018 
1019  /// Return face index of the e-th element adjacent to boundary b in region r
1020  int face_index_at_boundary_in_region(const unsigned& b,
1021  const unsigned& r,
1022  const unsigned& e) const
1023  {
1024  // Use a constant iterator here to keep function "const" overall
1025  std::map<unsigned, Vector<int>>::const_iterator it =
1026  Face_index_region_at_boundary[b].find(r);
1027  if (it != Face_index_region_at_boundary[b].end())
1028  {
1029  return (it->second)[e];
1030  }
1031  else
1032  {
1033  std::ostringstream error_message;
1034  error_message << "Face indices not set up for boundary " << b
1035  << " in region " << r << "\n";
1036  error_message << "This probably means that the boundary is not "
1037  "adjacent to region\n";
1038  throw OomphLibError(error_message.str(),
1039  OOMPH_CURRENT_FUNCTION,
1040  OOMPH_EXCEPTION_LOCATION);
1041  }
1042  }
1043 
1044 
1045  /// Return the number of regions specified by attributes
1046  unsigned nregion()
1047  {
1048  return Region_element_pt.size();
1049  }
1050 
1051  /// Return the number of elements in region r
1052  unsigned nregion_element(const unsigned& r)
1053  {
1054  unsigned entry = 0;
1055  bool found = false;
1056  unsigned n = Region_attribute.size();
1057  for (unsigned i = 0; i < n; i++)
1058  {
1059 #ifdef PARANOID
1060  double diff =
1061  fabs(Region_attribute[i] -
1062  static_cast<double>(static_cast<unsigned>(Region_attribute[i])));
1063  if (diff > 0.0)
1064  {
1065  std::ostringstream error_message;
1066  error_message << "Region attributes should be unsigneds because we \n"
1067  << "only use them to set region ids\n";
1068  throw OomphLibError(error_message.str(),
1069  OOMPH_CURRENT_FUNCTION,
1070  OOMPH_EXCEPTION_LOCATION);
1071  }
1072 #endif
1073  if (static_cast<unsigned>(Region_attribute[i]) == r)
1074  {
1075  entry = i;
1076  found = true;
1077  break;
1078  }
1079  }
1080  if (found)
1081  {
1082  return Region_element_pt[entry].size();
1083  }
1084  else
1085  {
1086  return 0;
1087  }
1088  }
1089 
1090  /// Return the i-th region attribute (here only used as the
1091  /// (assumed to be unsigned) region id
1092  double region_attribute(const unsigned& i)
1093  {
1094  return Region_attribute[i];
1095  }
1096 
1097  /// Return the e-th element in the r-th region
1098  FiniteElement* region_element_pt(const unsigned& r, const unsigned& e)
1099  {
1100  unsigned entry = 0;
1101  bool found = false;
1102  unsigned n = Region_attribute.size();
1103  for (unsigned i = 0; i < n; i++)
1104  {
1105 #ifdef PARANOID
1106  double diff =
1107  fabs(Region_attribute[i] -
1108  static_cast<double>(static_cast<unsigned>(Region_attribute[i])));
1109  if (diff > 0.0)
1110  {
1111  std::ostringstream error_message;
1112  error_message << "Region attributes should be unsigneds because we \n"
1113  << "only use the to set region ids\n";
1114  throw OomphLibError(error_message.str(),
1115  OOMPH_CURRENT_FUNCTION,
1116  OOMPH_EXCEPTION_LOCATION);
1117  }
1118 #endif
1119  if (static_cast<unsigned>(Region_attribute[i]) == r)
1120  {
1121  entry = i;
1122  found = true;
1123  break;
1124  }
1125  }
1126  if (found)
1127  {
1128  return Region_element_pt[entry][e];
1129  }
1130  else
1131  {
1132  return 0;
1133  }
1134  }
1135 
1136  /// Snap boundaries specified by the IDs listed in boundary_id to
1137  /// a quadratric surface, specified in the file
1138  /// quadratic_surface_file_name. This is usually used with vmtk-based
1139  /// meshes for which oomph-lib's xda to poly conversion code produces the
1140  /// files "quadratic_fsi_boundary.dat" and
1141  /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
1142  /// boundary (for the fluid and the solid) and the quadratic representation
1143  /// of the outer boundary of the solid. When used with these files, the flag
1144  /// switch_normal should be set to true when calling the function for the
1145  /// outer boundary of the solid. The DocInfo object can be used to label
1146  /// optional output files. (Uses directory and label).
1147  template<class ELEMENT>
1149  const Vector<unsigned>& boundary_id,
1150  const std::string& quadratic_surface_file_name,
1151  const bool& switch_normal,
1152  DocInfo& doc_info);
1153 
1154  /// Snap boundaries specified by the IDs listed in boundary_id to
1155  /// a quadratric surface, specified in the file
1156  /// quadratic_surface_file_name. This is usually used with vmtk-based
1157  /// meshes for which oomph-lib's xda to poly conversion code produces the
1158  /// files "quadratic_fsi_boundary.dat" and
1159  /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
1160  /// boundary (for the fluid and the solid) and the quadratic representation
1161  /// of the outer boundary of the solid. When used with these files, the flag
1162  /// switch_normal should be set to true when calling the function for the
1163  /// outer boundary of the solid.
1164  template<class ELEMENT>
1166  const Vector<unsigned>& boundary_id,
1167  const std::string& quadratic_surface_file_name,
1168  const bool& switch_normal)
1169  {
1170  // Dummy doc info
1171  DocInfo doc_info;
1172  doc_info.disable_doc();
1173  snap_to_quadratic_surface<ELEMENT>(
1174  boundary_id, quadratic_surface_file_name, switch_normal, doc_info);
1175  }
1176 
1177  /// Move the nodes on boundaries with associated GeomObjects so
1178  /// that they exactly coincide with the geometric object. This requires
1179  /// that the boundary coordinates are set up consistently.
1181 
1182 
1183  /// Non-Delaunay split elements that have three faces on a boundary
1184  /// into sons. Timestepper species timestepper for new nodes; defaults
1185  /// to to steady timestepper.
1186  template<class ELEMENT>
1188  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper);
1189 
1190 
1191  /// Setup lookup schemes which establish which elements are located
1192  /// next to mesh's boundaries (wrapper to suppress doc).
1194  {
1195  std::ofstream outfile;
1196  this->setup_boundary_element_info(outfile);
1197  }
1198 
1199 
1200  /// Setup lookup schemes which establish which elements are located
1201  /// next to mesh's boundaries. Doc in outfile (if it's open).
1202  void setup_boundary_element_info(std::ostream& outfile);
1203 
1204 
1205  protected:
1206  /// Vectors of vectors of elements in each region (note: this just
1207  /// stores them; the region IDs are contained in Region_attribute!)
1209 
1210  /// Vector of attributes associated with the elements in each region
1211  /// NOTE: double is enforced on us by tetgen. We use it as an unsigned
1212  /// to indicate the actual (zero-based) region ID
1214 
1215  /// Storage for elements adjacent to a boundary in a particular region
1218 
1219  /// Storage for the face index adjacent to a boundary in a particular
1220  /// region
1222 
1223  /// Faceted surface that defines outer boundaries
1225 
1226  /// Vector to faceted surfaces that define internal boundaries
1228 
1229  /// Reverse lookup scheme: Pointer to faceted surface (if any!)
1230  /// associated with boundary b
1231  std::map<unsigned, TetMeshFacetedSurface*> Tet_mesh_faceted_surface_pt;
1232 
1233  /// Reverse lookup scheme: Pointer to facet (if any!)
1234  /// associated with boundary b
1235  std::map<unsigned, TetMeshFacet*> Tet_mesh_facet_pt;
1236 
1237  /// Boundary coordinates of vertices in triangular facets
1238  /// associated with given boundary. Is only set up for triangular facets!
1239  std::map<unsigned, Vector<Vector<double>>>
1241 
1242  /// Timestepper used to build nodes
1244  };
1245 
1246 
1247  /// ///////////////////////////////////////////////////////////////////////////
1248  /// ///////////////////////////////////////////////////////////////////////////
1249  /// ///////////////////////////////////////////////////////////////////////////
1250 
1251 
1252  //###########################################################################
1253  // Templated member functions
1254  //###########################################################################
1255 
1256 
1257  //======================================================================
1258  /// Setup boundary coordinate on boundary b which is
1259  /// assumed to be planar. Boundary coordinates are the
1260  /// x-y coordinates in the plane of that boundary, with the
1261  /// x-axis along the line from the (lexicographically)
1262  /// "lower left" to the "upper right" node. The y axis
1263  /// is obtained by taking the cross-product of the positive
1264  /// x direction with the outer unit normal computed by
1265  /// the face elements (or its negative if switch_normal is set
1266  /// to true). Doc faces in output file (if it's open).
1267  ///
1268  /// Note 1: Setup of boundary coordinates is not done if the boundary in
1269  /// question turns out to be nonplanar.
1270  ///
1271  /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
1272  /// also
1273  /// store the boundary coordinates of its vertices. They are needed
1274  /// to interpolated intrinsic coordinates of an associated GeomObject
1275  /// (if any) into the interior.
1276  //======================================================================
1277  template<class ELEMENT>
1279  const bool& switch_normal,
1280  std::ofstream& outfile)
1281  {
1282  Node* lower_left_node_pt = 0;
1283  Node* upper_right_node_pt = 0;
1284 
1285  // Unit vector connecting lower left and upper right nodes
1286  Vector<double> b0(3);
1287 
1288  // ...and a vector normal to it
1289  Vector<double> b1(3);
1290 
1291 
1292  // Facet?
1293  TetMeshFacet* f_pt = 0;
1294  std::map<unsigned, TetMeshFacet*>::iterator it = Tet_mesh_facet_pt.find(b);
1295  if (it != Tet_mesh_facet_pt.end())
1296  {
1297  f_pt = (*it).second;
1298  }
1299 
1300  // std::cout << "Debug " << b; f_pt->output(std::cout);
1301 
1302  // Number of vertices
1303  unsigned nv = 0;
1304  if (f_pt != 0)
1305  {
1306  nv = f_pt->nvertex();
1307  }
1308 
1309  // Check for planarity and bail out if nodes are not in the same plane
1310  bool vertices_are_in_same_plane = true;
1311  for (unsigned do_for_real = 0; do_for_real < 2; do_for_real++)
1312  {
1313  // Temporary storage for face elements
1314  Vector<FiniteElement*> face_el_pt;
1315 
1316  // Loop over all elements on boundaries
1317  unsigned nel = this->nboundary_element(b);
1318 
1319  if (nel > 0)
1320  {
1321  // Loop over the bulk elements adjacent to boundary b
1322  for (unsigned e = 0; e < nel; e++)
1323  {
1324  // Get pointer to the bulk element that is adjacent to boundary b
1325  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
1326 
1327  // Find the index of the face of element e along boundary b
1328  int face_index = this->face_index_at_boundary(b, e);
1329 
1330  // Create new face element
1331  face_el_pt.push_back(
1332  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index));
1333 
1334  // Output faces?
1335  if (outfile.is_open())
1336  {
1337  face_el_pt[face_el_pt.size() - 1]->output(outfile);
1338  }
1339  }
1340 
1341  // Loop over all nodes to find the lower left and upper right ones
1342  lower_left_node_pt = this->boundary_node_pt(b, 0);
1343  upper_right_node_pt = this->boundary_node_pt(b, 0);
1344 
1345  // Loop over all nodes on the boundary
1346  unsigned nnod = this->nboundary_node(b);
1347  for (unsigned j = 0; j < nnod; j++)
1348  {
1349  // Get node
1350  Node* nod_pt = this->boundary_node_pt(b, j);
1351 
1352  // Primary criterion for lower left: Does it have a lower
1353  // z-coordinate?
1354  if (nod_pt->x(2) < lower_left_node_pt->x(2))
1355  {
1356  lower_left_node_pt = nod_pt;
1357  }
1358  // ...or is it a draw?
1359  else if (nod_pt->x(2) == lower_left_node_pt->x(2))
1360  {
1361  // If it's a draw: Does it have a lower y-coordinate?
1362  if (nod_pt->x(1) < lower_left_node_pt->x(1))
1363  {
1364  lower_left_node_pt = nod_pt;
1365  }
1366  // ...or is it a draw?
1367  else if (nod_pt->x(1) == lower_left_node_pt->x(1))
1368  {
1369  // If it's a draw: Does it have a lower x-coordinate?
1370  if (nod_pt->x(0) < lower_left_node_pt->x(0))
1371  {
1372  lower_left_node_pt = nod_pt;
1373  }
1374  }
1375  }
1376 
1377  // Primary criterion for upper right: Does it have a higher
1378  // z-coordinate?
1379  if (nod_pt->x(2) > upper_right_node_pt->x(2))
1380  {
1381  upper_right_node_pt = nod_pt;
1382  }
1383  // ...or is it a draw?
1384  else if (nod_pt->x(2) == upper_right_node_pt->x(2))
1385  {
1386  // If it's a draw: Does it have a higher y-coordinate?
1387  if (nod_pt->x(1) > upper_right_node_pt->x(1))
1388  {
1389  upper_right_node_pt = nod_pt;
1390  }
1391  // ...or is it a draw?
1392  else if (nod_pt->x(1) == upper_right_node_pt->x(1))
1393  {
1394  // If it's a draw: Does it have a higher x-coordinate?
1395  if (nod_pt->x(0) > upper_right_node_pt->x(0))
1396  {
1397  upper_right_node_pt = nod_pt;
1398  }
1399  }
1400  }
1401  }
1402 
1403  // Prepare storage for boundary coordinates
1404  Vector<double> zeta(2);
1405  /*std::cout << upper_right_node_pt->x(0) << " "
1406  << upper_right_node_pt->x(1) << " "
1407  << upper_right_node_pt->x(2) << " L ";
1408  std::cout << lower_left_node_pt->x(0) << " "
1409  << lower_left_node_pt->x(1) << " "
1410  << lower_left_node_pt->x(2) << "\n ";*/
1411 
1412 
1413  // Unit vector connecting lower left and upper right nodes
1414  b0[0] = upper_right_node_pt->x(0) - lower_left_node_pt->x(0);
1415  b0[1] = upper_right_node_pt->x(1) - lower_left_node_pt->x(1);
1416  b0[2] = upper_right_node_pt->x(2) - lower_left_node_pt->x(2);
1417 
1418  // Normalise
1419  double inv_length =
1420  1.0 / sqrt(b0[0] * b0[0] + b0[1] * b0[1] + b0[2] * b0[2]);
1421  b0[0] *= inv_length;
1422  b0[1] *= inv_length;
1423  b0[2] *= inv_length;
1424 
1425  // std::cout << "B0 ";
1426  // std::cout << b0[0] << " " << b0[1] << " " << b0[2] << "\n";
1427 
1428  // Get (outer) unit normal to first face element
1429  Vector<double> normal(3);
1430  Vector<double> s(2, 0.0);
1431  if (nv != 3)
1432  {
1433  dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[0])
1434  ->outer_unit_normal(s, normal);
1435  }
1436  else
1437  {
1438  double t1x = f_pt->vertex_pt(1)->x(0) - f_pt->vertex_pt(0)->x(0);
1439 
1440  double t1y = f_pt->vertex_pt(1)->x(1) - f_pt->vertex_pt(0)->x(1);
1441 
1442  double t1z = f_pt->vertex_pt(1)->x(2) - f_pt->vertex_pt(0)->x(2);
1443 
1444  double t2x = f_pt->vertex_pt(2)->x(0) - f_pt->vertex_pt(0)->x(0);
1445 
1446  double t2y = f_pt->vertex_pt(2)->x(1) - f_pt->vertex_pt(0)->x(1);
1447 
1448  double t2z = f_pt->vertex_pt(2)->x(2) - f_pt->vertex_pt(0)->x(2);
1449 
1450  normal[0] = t1y * t2z - t1z * t2y;
1451  normal[1] = t1z * t2x - t1x * t2z;
1452  normal[2] = t1x * t2y - t1y * t2x;
1453  double inv_length =
1454  1.0 / sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
1455  normal[2] * normal[2]);
1456  normal[0] *= inv_length;
1457  normal[1] *= inv_length;
1458  normal[2] *= inv_length;
1459  }
1460 
1461  if (switch_normal)
1462  {
1463  normal[0] = -normal[0];
1464  normal[1] = -normal[1];
1465  normal[2] = -normal[2];
1466  }
1467 
1468  // std::cout << "Normal ";
1469  // std::cout << normal[0] << " " << normal[1] << " " << normal[2] <<
1470  // "\n";
1471 
1472 
1473  // Check that all elements have the same normal
1474  for (unsigned e = 0; e < nel; e++)
1475  {
1476  // Get (outer) unit normal to face element
1477  Vector<double> my_normal(3);
1478  dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[e])
1479  ->outer_unit_normal(s, my_normal);
1480 
1481  // Dot product should be one!
1482  double dot_prod = normal[0] * my_normal[0] +
1483  normal[1] * my_normal[1] + normal[2] * my_normal[2];
1484 
1485 
1486  double error = abs(dot_prod) - 1.0;
1487  if (abs(error) > Tolerance_for_boundary_finding)
1488  {
1489  vertices_are_in_same_plane = false;
1490  }
1491  }
1492 
1493  // Bail out?
1494  if (vertices_are_in_same_plane)
1495  {
1496  // Cross-product to get second in-plane vector, normal to b0
1497  b1[0] = b0[1] * normal[2] - b0[2] * normal[1];
1498  b1[1] = b0[2] * normal[0] - b0[0] * normal[2];
1499  b1[2] = b0[0] * normal[1] - b0[1] * normal[0];
1500 
1501  // std::cout << "B1 ";
1502  // std::cout << b1[0] << " " << b1[1] << " " << b1[2] << "\n";
1503 
1504 
1505  // Assign boundary coordinates: projection onto the axes
1506  for (unsigned j = 0; j < nnod; j++)
1507  {
1508  // Get node
1509  Node* nod_pt = this->boundary_node_pt(b, j);
1510 
1511  // Difference vector to lower left corner
1512  Vector<double> delta(3);
1513  delta[0] = nod_pt->x(0) - lower_left_node_pt->x(0);
1514  delta[1] = nod_pt->x(1) - lower_left_node_pt->x(1);
1515  delta[2] = nod_pt->x(2) - lower_left_node_pt->x(2);
1516 
1517  /*std::cout << j << ": "
1518  << nod_pt->x(0) << " " << nod_pt->x(1)
1519  << " " << nod_pt->x(2) << " Delta ";
1520  std::cout << delta[0] << " " << delta[1] << " " << delta[2] << "\n";
1521  */
1522 
1523  // Project
1524  zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1525  zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1526 
1527  // Check:
1528  double error = pow(lower_left_node_pt->x(0) + zeta[0] * b0[0] +
1529  zeta[1] * b1[0] - nod_pt->x(0),
1530  2) +
1531  pow(lower_left_node_pt->x(1) + zeta[0] * b0[1] +
1532  zeta[1] * b1[1] - nod_pt->x(1),
1533  2) +
1534  pow(lower_left_node_pt->x(2) + zeta[0] * b0[2] +
1535  zeta[1] * b1[2] - nod_pt->x(2),
1536  2);
1537 
1538  if (sqrt(error) > Tolerance_for_boundary_finding)
1539  {
1540  std::ostringstream error_message;
1541  error_message
1542  << "Error in projection of boundary coordinates = "
1543  << sqrt(error) << " > Tolerance_for_boundary_finding = "
1544  << Tolerance_for_boundary_finding << std::endl
1545  << "nv = " << nv << std::endl
1546  << "Lower left: " << lower_left_node_pt->x(0) << " "
1547  << lower_left_node_pt->x(1) << " " << lower_left_node_pt->x(2)
1548  << " " << std::endl
1549  << "Upper right: " << upper_right_node_pt->x(0) << " "
1550  << upper_right_node_pt->x(1) << " " << upper_right_node_pt->x(2)
1551  << " " << std::endl
1552  << "Nodes: ";
1553  for (unsigned j = 0; j < nnod; j++)
1554  {
1555  // Get node
1556  Node* nod_pt = this->boundary_node_pt(b, j);
1557  error_message << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1558  << nod_pt->x(2) << " " << std::endl;
1559  }
1560  throw OomphLibError(error_message.str(),
1561  OOMPH_CURRENT_FUNCTION,
1562  OOMPH_EXCEPTION_LOCATION);
1563  }
1564 
1565  // Set boundary coordinate
1566  if (do_for_real == 1)
1567  {
1568  nod_pt->set_coordinates_on_boundary(b, zeta);
1569  }
1570  }
1571  } // End of if vertices are in the same plane
1572  }
1573 
1574 
1575  // Indicate that boundary coordinate has been set up
1576  if (do_for_real == 1)
1577  {
1578  Boundary_coordinate_exists[b] = true;
1579 
1580  // Facet associated with this boundary?
1581  if (f_pt != 0)
1582  {
1583  // Triangular facet: Set coordinates at vertices
1584  if (nv == 3)
1585  {
1587  for (unsigned j = 0; j < 3; j++)
1588  {
1589  // Two surface coordinates
1591 
1592  // Difference vector to lower left corner
1593  Vector<double> delta(3);
1594  delta[0] = f_pt->vertex_pt(j)->x(0) - lower_left_node_pt->x(0);
1595  delta[1] = f_pt->vertex_pt(j)->x(1) - lower_left_node_pt->x(1);
1596  delta[2] = f_pt->vertex_pt(j)->x(2) - lower_left_node_pt->x(2);
1597 
1598  // Project
1599  Vector<double> zeta(2);
1600  zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1601  zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1602 
1603  for (unsigned ii = 0; ii < 2; ii++)
1604  {
1606  zeta[ii];
1607  }
1608  }
1609  }
1610  }
1611  }
1612 
1613  // Cleanup
1614  unsigned n = face_el_pt.size();
1615  for (unsigned e = 0; e < n; e++)
1616  {
1617  delete face_el_pt[e];
1618  }
1619 
1620  // Bail out?
1621  if (!vertices_are_in_same_plane)
1622  {
1623  /* oomph_info << "Vertices on boundary " << b */
1624  /* << " are not in the same plane; bailing out\n"; */
1625  return;
1626  }
1627  }
1628  }
1629 
1630 
1631  //======================================================================
1632  /// Snap boundaries specified by the IDs listed in boundary_id to
1633  /// a quadratric surface, specified in the file
1634  /// quadratic_surface_file_name. This is usually used with vmtk-based
1635  /// meshes for which oomph-lib's xda to poly conversion code produces the
1636  /// files "quadratic_fsi_boundary.dat" and
1637  /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
1638  /// boundary (for the fluid and the solid) and the quadratic representation of
1639  /// the outer boundary of the solid. When used with these files, the flag
1640  /// switch_normal should be set to true when calling the function for the
1641  /// outer boundary of the solid. The DocInfo object can be used to label
1642  /// optional output files. (Uses directory and label).
1643  //======================================================================
1644  template<class ELEMENT>
1646  const Vector<unsigned>& boundary_id,
1647  const std::string& quadratic_surface_file_name,
1648  const bool& switch_normal,
1649  DocInfo& doc_info)
1650  {
1651  // Aux storage for processing input
1652  char dummy[101];
1653 
1654  // Prepare to doc nodes that couldn't be snapped
1655  std::ofstream the_file_non_snapped;
1656  std::string non_snapped_filename =
1657  "non_snapped_nodes_" + doc_info.label() + ".dat";
1658 
1659  // Read the number of nodes and elements (quadratic facets)
1660  std::ifstream infile(quadratic_surface_file_name.c_str(),
1661  std::ios_base::in);
1662  unsigned n_node;
1663  infile >> n_node;
1664 
1665  // Ignore rest of line
1666  infile.getline(dummy, 100);
1667 
1668  // Number of quadratic facets
1669  unsigned nel;
1670  infile >> nel;
1671 
1672  // Ignore rest of line
1673  infile.getline(dummy, 100);
1674 
1675  // Check that the number of elements matches
1676  if (nel != boundary_id.size())
1677  {
1678  std::ostringstream error_message;
1679  error_message << "Number of quadratic facets specified in "
1680  << quadratic_surface_file_name << "is: " << nel
1681  << "\nThis doesn't match the number of planar boundaries \n"
1682  << "specified in boundary_id which is "
1683  << boundary_id.size() << std::endl;
1684  throw OomphLibError(
1685  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1686  }
1687 
1688  // Temporary storage for face elements
1690 
1691  // Loop over quadratic face elements -- one for each facet
1692  for (unsigned e = 0; e < nel; e++)
1693  {
1694  face_el_pt.push_back(new FreeStandingFaceElement<ELEMENT>);
1695  }
1696 
1697 
1698  // Now build nodes
1699  unsigned n_dim = 3;
1700  unsigned n_position_type = 1;
1701  unsigned initial_n_value = 0;
1702  Vector<Node*> nod_pt(n_node);
1703  unsigned node_nmbr;
1704  std::map<unsigned, unsigned> node_number;
1705  std::ofstream node_file;
1706  for (unsigned j = 0; j < n_node; j++)
1707  {
1708  nod_pt[j] =
1709  new BoundaryNode<Node>(n_dim, n_position_type, initial_n_value);
1710  infile >> nod_pt[j]->x(0);
1711  infile >> nod_pt[j]->x(1);
1712  infile >> nod_pt[j]->x(2);
1713  infile >> node_nmbr;
1714  node_number[node_nmbr] = j;
1715  }
1716 
1717 
1718  // Now assign nodes to elements -- each element represents
1719  // distinct boundary; assign enumeration as specified by
1720  // boundary_id.
1721  for (unsigned e = 0; e < nel; e++)
1722  {
1723  unsigned nnod_el = face_el_pt[e]->nnode();
1724  unsigned j_global;
1725  for (unsigned j = 0; j < nnod_el; j++)
1726  {
1727  infile >> j_global;
1728  face_el_pt[e]->node_pt(j) = nod_pt[node_number[j_global]];
1729  face_el_pt[e]->node_pt(j)->add_to_boundary(boundary_id[e]);
1730  }
1731  face_el_pt[e]->set_boundary_number_in_bulk_mesh(boundary_id[e]);
1732  face_el_pt[e]->set_nodal_dimension(3);
1733  }
1734 
1735 
1736  // Setup boundary coordinates for each facet, using
1737  // the same strategy as for the simplex boundaries
1738  // (there's some code duplication here but it doesn't
1739  // seem worth it to rationlise this as the interface would
1740  // be pretty clunky).
1741  for (unsigned e = 0; e < nel; e++)
1742  {
1743  Vector<Vector<double>> vertex_boundary_coord(3);
1744 
1745  // Loop over all nodes to find the lower left and upper right ones
1746  Node* lower_left_node_pt = face_el_pt[e]->node_pt(0);
1747  Node* upper_right_node_pt = face_el_pt[e]->node_pt(0);
1748 
1749  // Loop over all vertex nodes
1750  for (unsigned j = 0; j < 3; j++)
1751  {
1752  // Get node
1753  Node* nod_pt = face_el_pt[e]->node_pt(j);
1754 
1755  // Primary criterion for lower left: Does it have a lower z-coordinate?
1756  if (nod_pt->x(2) < lower_left_node_pt->x(2))
1757  {
1758  lower_left_node_pt = nod_pt;
1759  }
1760  // ...or is it a draw?
1761  else if (nod_pt->x(2) == lower_left_node_pt->x(2))
1762  {
1763  // If it's a draw: Does it have a lower y-coordinate?
1764  if (nod_pt->x(1) < lower_left_node_pt->x(1))
1765  {
1766  lower_left_node_pt = nod_pt;
1767  }
1768  // ...or is it a draw?
1769  else if (nod_pt->x(1) == lower_left_node_pt->x(1))
1770  {
1771  // If it's a draw: Does it have a lower x-coordinate?
1772  if (nod_pt->x(0) < lower_left_node_pt->x(0))
1773  {
1774  lower_left_node_pt = nod_pt;
1775  }
1776  }
1777  }
1778 
1779  // Primary criterion for upper right: Does it have a higher
1780  // z-coordinate?
1781  if (nod_pt->x(2) > upper_right_node_pt->x(2))
1782  {
1783  upper_right_node_pt = nod_pt;
1784  }
1785  // ...or is it a draw?
1786  else if (nod_pt->x(2) == upper_right_node_pt->x(2))
1787  {
1788  // If it's a draw: Does it have a higher y-coordinate?
1789  if (nod_pt->x(1) > upper_right_node_pt->x(1))
1790  {
1791  upper_right_node_pt = nod_pt;
1792  }
1793  // ...or is it a draw?
1794  else if (nod_pt->x(1) == upper_right_node_pt->x(1))
1795  {
1796  // If it's a draw: Does it have a higher x-coordinate?
1797  if (nod_pt->x(0) > upper_right_node_pt->x(0))
1798  {
1799  upper_right_node_pt = nod_pt;
1800  }
1801  }
1802  }
1803  }
1804 
1805  // Prepare storage for boundary coordinates
1806  Vector<double> zeta(2);
1807 
1808  // Unit vector connecting lower left and upper right nodes
1809  Vector<double> b0(3);
1810  b0[0] = upper_right_node_pt->x(0) - lower_left_node_pt->x(0);
1811  b0[1] = upper_right_node_pt->x(1) - lower_left_node_pt->x(1);
1812  b0[2] = upper_right_node_pt->x(2) - lower_left_node_pt->x(2);
1813 
1814  // Normalise
1815  double inv_length =
1816  1.0 / sqrt(b0[0] * b0[0] + b0[1] * b0[1] + b0[2] * b0[2]);
1817  b0[0] *= inv_length;
1818  b0[1] *= inv_length;
1819  b0[2] *= inv_length;
1820 
1821  // Get (outer) unit normal to face element -- note that
1822  // with the enumeration chosen in oomph-lib's xda to poly
1823  // conversion code the sign below is correct for the outer
1824  // unit normal on the FSI interface.
1825  Vector<double> tang1(3);
1826  Vector<double> tang2(3);
1827  Vector<double> normal(3);
1828  tang1[0] =
1829  face_el_pt[e]->node_pt(1)->x(0) - face_el_pt[e]->node_pt(0)->x(0);
1830  tang1[1] =
1831  face_el_pt[e]->node_pt(1)->x(1) - face_el_pt[e]->node_pt(0)->x(1);
1832  tang1[2] =
1833  face_el_pt[e]->node_pt(1)->x(2) - face_el_pt[e]->node_pt(0)->x(2);
1834  tang2[0] =
1835  face_el_pt[e]->node_pt(2)->x(0) - face_el_pt[e]->node_pt(0)->x(0);
1836  tang2[1] =
1837  face_el_pt[e]->node_pt(2)->x(1) - face_el_pt[e]->node_pt(0)->x(1);
1838  tang2[2] =
1839  face_el_pt[e]->node_pt(2)->x(2) - face_el_pt[e]->node_pt(0)->x(2);
1840  normal[0] = -(tang1[1] * tang2[2] - tang1[2] * tang2[1]);
1841  normal[1] = -(tang1[2] * tang2[0] - tang1[0] * tang2[2]);
1842  normal[2] = -(tang1[0] * tang2[1] - tang1[1] * tang2[0]);
1843 
1844  // Normalise
1845  inv_length = 1.0 / sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
1846  normal[2] * normal[2]);
1847  normal[0] *= inv_length;
1848  normal[1] *= inv_length;
1849  normal[2] *= inv_length;
1850 
1851  // Change direction -- usually for outer boundary of solid
1852  if (switch_normal)
1853  {
1854  normal[0] = -normal[0];
1855  normal[1] = -normal[1];
1856  normal[2] = -normal[2];
1857  }
1858 
1859  // Cross-product to get second in-plane vector, normal to b0
1860  Vector<double> b1(3);
1861  b1[0] = b0[1] * normal[2] - b0[2] * normal[1];
1862  b1[1] = b0[2] * normal[0] - b0[0] * normal[2];
1863  b1[2] = b0[0] * normal[1] - b0[1] * normal[0];
1864 
1865  // Assign boundary coordinates for vertex nodes: projection onto the axes
1866  for (unsigned j = 0; j < 3; j++)
1867  {
1868  // Get node
1869  Node* nod_pt = face_el_pt[e]->node_pt(j);
1870 
1871  // Difference vector to lower left corner
1872  Vector<double> delta(3);
1873  delta[0] = nod_pt->x(0) - lower_left_node_pt->x(0);
1874  delta[1] = nod_pt->x(1) - lower_left_node_pt->x(1);
1875  delta[2] = nod_pt->x(2) - lower_left_node_pt->x(2);
1876 
1877  // Project
1878  zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1879  zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1880 
1881  vertex_boundary_coord[j].resize(2);
1882  vertex_boundary_coord[j][0] = zeta[0];
1883  vertex_boundary_coord[j][1] = zeta[1];
1884 
1885 #ifdef PARANOID
1886 
1887  // Check:
1888  double error = pow(lower_left_node_pt->x(0) + zeta[0] * b0[0] +
1889  zeta[1] * b1[0] - nod_pt->x(0),
1890  2) +
1891  pow(lower_left_node_pt->x(1) + zeta[0] * b0[1] +
1892  zeta[1] * b1[1] - nod_pt->x(1),
1893  2) +
1894  pow(lower_left_node_pt->x(2) + zeta[0] * b0[2] +
1895  zeta[1] * b1[2] - nod_pt->x(2),
1896  2);
1897 
1898  if (sqrt(error) > Tolerance_for_boundary_finding)
1899  {
1900  std::ostringstream error_message;
1901  error_message
1902  << "Error in setup of boundary coordinate " << sqrt(error)
1903  << std::endl
1904  << "exceeds tolerance specified by the static member data\n"
1905  << "TetMeshBase::Tolerance_for_boundary_finding = "
1906  << Tolerance_for_boundary_finding << std::endl
1907  << "This usually means that the boundary is not planar.\n\n"
1908  << "You can suppress this error message by recompiling \n"
1909  << "recompiling without PARANOID or by changing the tolerance.\n";
1910  throw OomphLibError(error_message.str(),
1911  OOMPH_CURRENT_FUNCTION,
1912  OOMPH_EXCEPTION_LOCATION);
1913  }
1914 #endif
1915 
1916  // Set boundary coordinate
1917  nod_pt->set_coordinates_on_boundary(boundary_id[e], zeta);
1918  }
1919 
1920  // Assign boundary coordinates of three midside nodes by linear
1921  // interpolation (corresponding to a flat facet)
1922 
1923  // Node 3 is between 0 and 1
1924  zeta[0] =
1925  0.5 * (vertex_boundary_coord[0][0] + vertex_boundary_coord[1][0]);
1926  zeta[1] =
1927  0.5 * (vertex_boundary_coord[0][1] + vertex_boundary_coord[1][1]);
1928  face_el_pt[e]->node_pt(3)->set_coordinates_on_boundary(boundary_id[e],
1929  zeta);
1930 
1931  // Node 4 is between 1 and 2
1932  zeta[0] =
1933  0.5 * (vertex_boundary_coord[1][0] + vertex_boundary_coord[2][0]);
1934  zeta[1] =
1935  0.5 * (vertex_boundary_coord[1][1] + vertex_boundary_coord[2][1]);
1936  face_el_pt[e]->node_pt(4)->set_coordinates_on_boundary(boundary_id[e],
1937  zeta);
1938 
1939  // Node 5 is between 2 and 0
1940  zeta[0] =
1941  0.5 * (vertex_boundary_coord[2][0] + vertex_boundary_coord[0][0]);
1942  zeta[1] =
1943  0.5 * (vertex_boundary_coord[2][1] + vertex_boundary_coord[0][1]);
1944  face_el_pt[e]->node_pt(5)->set_coordinates_on_boundary(boundary_id[e],
1945  zeta);
1946  }
1947 
1948 
1949  // Loop over elements/facets = boundaries to snap
1950  bool success = true;
1951  for (unsigned b = 0; b < nel; b++)
1952  {
1953  // Doc boundary coordinates on quadratic patches
1954  std::ofstream the_file;
1955  std::ofstream the_file_before;
1956  std::ofstream the_file_after;
1957  if (doc_info.is_doc_enabled())
1958  {
1959  std::ostringstream filename;
1960  filename << doc_info.directory() << "/quadratic_coordinates_"
1961  << doc_info.label() << b << ".dat";
1962  the_file.open(filename.str().c_str());
1963 
1964  std::ostringstream filename1;
1965  filename1 << doc_info.directory() << "/quadratic_nodes_before_"
1966  << doc_info.label() << b << ".dat";
1967  the_file_before.open(filename1.str().c_str());
1968 
1969  std::ostringstream filename2;
1970  filename2 << doc_info.directory() << "/quadratic_nodes_after_"
1971  << doc_info.label() << b << ".dat";
1972  the_file_after.open(filename2.str().c_str());
1973  }
1974 
1975  // Cast the element pointer
1976  FreeStandingFaceElement<ELEMENT>* el_pt = face_el_pt[b];
1977 
1978  // Doc boundary coordinate on quadratic facet representation
1979  if (doc_info.is_doc_enabled())
1980  {
1981  Vector<double> s(2);
1982  Vector<double> zeta(2);
1983  Vector<double> x(3);
1984  unsigned n_plot = 3;
1985  the_file << el_pt->tecplot_zone_string(n_plot);
1986 
1987  // Loop over plot points
1988  unsigned num_plot_points = el_pt->nplot_points(n_plot);
1989  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1990  {
1991  // Get local coordinates of plot point
1992  el_pt->get_s_plot(iplot, n_plot, s);
1993  el_pt->interpolated_zeta(s, zeta);
1994  el_pt->interpolated_x(s, x);
1995  for (unsigned i = 0; i < 3; i++)
1996  {
1997  the_file << x[i] << " ";
1998  }
1999  for (unsigned i = 0; i < 2; i++)
2000  {
2001  the_file << zeta[i] << " ";
2002  }
2003  the_file << std::endl;
2004  }
2005  el_pt->write_tecplot_zone_footer(the_file, n_plot);
2006 
2007  // std::cout << "Facet " << b << " corresponds to quadratic
2008  // boundary "
2009  // << boundary_id[b] << " which contains "
2010  // << this->nboundary_element(boundary_id[b])
2011  // << " element[s] " << std::endl;
2012  }
2013 
2014 
2015  // Loop over bulk elements that are adjacent to quadratic boundary
2016  Vector<double> boundary_zeta(2);
2017  Vector<double> quadratic_coordinates(3);
2018  GeomObject* geom_obj_pt = 0;
2019  Vector<double> s_geom_obj(2);
2020  unsigned nel = this->nboundary_element(boundary_id[b]);
2021  for (unsigned e = 0; e < nel; e++)
2022  {
2023  // Get pointer to the bulk element that is adjacent to boundary
2024  FiniteElement* bulk_elem_pt =
2025  this->boundary_element_pt(boundary_id[b], e);
2026 
2027  // Loop over nodes
2028  unsigned nnod = bulk_elem_pt->nnode();
2029  for (unsigned j = 0; j < nnod; j++)
2030  {
2031  Node* nod_pt = bulk_elem_pt->node_pt(j);
2032  if (nod_pt->is_on_boundary(boundary_id[b]))
2033  {
2034  nod_pt->get_coordinates_on_boundary(boundary_id[b], boundary_zeta);
2035 
2036  // Doc it?
2037  if (doc_info.is_doc_enabled())
2038  {
2039  the_file_before << nod_pt->x(0) << " " << nod_pt->x(1) << " "
2040  << nod_pt->x(2) << " " << boundary_zeta[0] << " "
2041  << boundary_zeta[1] << " " << std::endl;
2042  }
2043 
2044  // Find local coordinate in quadratic facet
2045  el_pt->locate_zeta(boundary_zeta, geom_obj_pt, s_geom_obj);
2046  if (geom_obj_pt != 0)
2047  {
2048  geom_obj_pt->position(s_geom_obj, quadratic_coordinates);
2049  nod_pt->x(0) = quadratic_coordinates[0];
2050  nod_pt->x(1) = quadratic_coordinates[1];
2051  nod_pt->x(2) = quadratic_coordinates[2];
2052  }
2053  else
2054  {
2055  // Get ready to bail out below...
2056  success = false;
2057 
2058  std::ostringstream error_message;
2059  error_message
2060  << "Warning: Couldn't find GeomObject during snapping to\n"
2061  << "quadratic surface for boundary " << boundary_id[b]
2062  << ". I'm leaving the node where it was. Will bail out "
2063  "later.\n";
2064  OomphLibWarning(error_message.str(),
2065  "TetgenMesh::snap_to_quadratic_surface()",
2066  OOMPH_EXCEPTION_LOCATION);
2067  if (!the_file_non_snapped.is_open())
2068  {
2069  the_file_non_snapped.open(non_snapped_filename.c_str());
2070  }
2071  the_file_non_snapped << nod_pt->x(0) << " " << nod_pt->x(1) << " "
2072  << nod_pt->x(2) << " " << boundary_zeta[0]
2073  << " " << boundary_zeta[1] << " "
2074  << std::endl;
2075  }
2076 
2077  // Doc it?
2078  if (doc_info.is_doc_enabled())
2079  {
2080  the_file_after << nod_pt->x(0) << " " << nod_pt->x(1) << " "
2081  << nod_pt->x(2) << " " << boundary_zeta[0] << " "
2082  << boundary_zeta[1] << " " << std::endl;
2083  }
2084  }
2085  }
2086  }
2087 
2088  // Close doc file
2089  the_file.close();
2090  the_file_before.close();
2091  the_file_after.close();
2092  }
2093 
2094  // Bail out?
2095  if (!success)
2096  {
2097  std::ostringstream error_message;
2098  error_message
2099  << "Warning: Couldn't find GeomObject during snapping to\n"
2100  << "quadratic surface. Bailing out.\n"
2101  << "Nodes that couldn't be snapped are contained in \n"
2102  << "file: " << non_snapped_filename << ".\n"
2103  << "This problem may arise because the switch_normal flag was \n"
2104  << "set wrongly.\n";
2105  throw OomphLibError(
2106  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2107  }
2108 
2109  // Close
2110  if (!the_file_non_snapped.is_open())
2111  {
2112  the_file_non_snapped.close();
2113  }
2114 
2115  // Kill auxiliary FreeStandingFaceElements
2116  for (unsigned e = 0; e < nel; e++)
2117  {
2118  delete face_el_pt[e];
2119  face_el_pt[e] = 0;
2120  }
2121 
2122  // Kill boundary nodes
2123  unsigned nn = nod_pt.size();
2124  for (unsigned j = 0; j < nn; j++)
2125  {
2126  delete nod_pt[j];
2127  }
2128  }
2129 
2130 
2131  //========================================================================
2132  /// Non-delaunay split elements that have three faces on a boundary
2133  /// into sons.
2134  //========================================================================
2135  template<class ELEMENT>
2137  {
2138  // Setup boundary lookup scheme if required
2140  {
2142  }
2143 
2144  // Find out how many nodes we have along each element edge
2145  unsigned nnode_1d = finite_element_pt(0)->nnode_1d();
2146  // Find out the total number of nodes
2147  unsigned nnode = this->finite_element_pt(0)->nnode();
2148 
2149  // At the moment we're only able to deal with nnode_1d=2 or 3.
2150  if ((nnode_1d != 2) && (nnode_1d != 3))
2151  {
2152  std::ostringstream error_message;
2153  error_message << "Mesh generation from tetgen currently only works\n";
2154  error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
2155  error_message << "for nnode_1d=" << nnode_1d << std::endl;
2156 
2157  throw OomphLibError(
2158  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2159  }
2160 
2161  // Map to store how many boundaries elements are located on
2162  std::map<FiniteElement*, unsigned> count;
2163 
2164  // Loop over all boundaries
2165  unsigned nb = this->nboundary();
2166  for (unsigned b = 0; b < nb; b++)
2167  {
2168  // Loop over elements next to that boundary
2169  unsigned nel = this->nboundary_element(b);
2170  for (unsigned e = 0; e < nel; e++)
2171  {
2172  // Get pointer to element
2173  FiniteElement* el_pt = boundary_element_pt(b, e);
2174 
2175  // Bump up counter
2176  count[el_pt]++;
2177  }
2178  }
2179 
2180  // To avoid having to check the map for all elements (which will
2181  // inflate it to the size of all elements!), move offending elements
2182  // into set
2183  std::set<FiniteElement*> elements_to_be_split;
2184  for (std::map<FiniteElement*, unsigned>::iterator it = count.begin();
2185  it != count.end();
2186  it++)
2187  {
2188  // Get pointer to element: Key
2189  FiniteElement* el_pt = it->first;
2190 
2191  // Does it have to be split?
2192  if (it->second > 2)
2193  {
2194  elements_to_be_split.insert(el_pt);
2195  }
2196  }
2197 
2198  // Vector for retained or newly built elements
2199  unsigned nel = this->nelement();
2200  Vector<FiniteElement*> new_or_retained_el_pt;
2201  new_or_retained_el_pt.reserve(nel);
2202 
2203  // Map which returns the 4 newly created elements for each old corner
2204  // element
2205  std::map<FiniteElement*, Vector<FiniteElement*>>
2206  old_to_new_corner_element_map;
2207 
2208  // Now loop over all elements
2209  for (unsigned e = 0; e < nel; e++)
2210  {
2211  // Get pointer to element
2212  FiniteElement* el_pt = this->finite_element_pt(e);
2213 
2214  // Does it have to be split? I.e. is it in the set?
2215  std::set<FiniteElement*>::iterator it = std::find(
2216  elements_to_be_split.begin(), elements_to_be_split.end(), el_pt);
2217 
2218  // It's not in the set, so iterator has reached end
2219  if (it == elements_to_be_split.end())
2220  {
2221  // Carry it across
2222  new_or_retained_el_pt.push_back(el_pt);
2223  }
2224  // It's in the set of elements to be split
2225  else
2226  {
2227  // Storage for new nodes -- Note: All new nodes are interior and
2228  // therefore cannot be boundary nodes!
2229  Node* node0_pt = 0;
2230  Node* node1_pt = 0;
2231  Node* node2_pt = 0;
2232  Node* node3_pt = 0;
2233  Node* node4_pt = 0;
2234  Node* node5_pt = 0;
2235  Node* node6_pt = 0;
2236  Node* node7_pt = 0;
2237  Node* node8_pt = 0;
2238  Node* node9_pt = 0;
2239  Node* node10_pt = 0;
2240 
2241  // Create first new element
2242  FiniteElement* el1_pt = new ELEMENT;
2243 
2244  // Copy existing nodes
2245  el1_pt->node_pt(0) = el_pt->node_pt(0);
2246  el1_pt->node_pt(1) = el_pt->node_pt(1);
2247  el1_pt->node_pt(3) = el_pt->node_pt(3);
2248  if (nnode_1d == 3)
2249  {
2250  el1_pt->node_pt(4) = el_pt->node_pt(4);
2251  el1_pt->node_pt(6) = el_pt->node_pt(6);
2252  el1_pt->node_pt(9) = el_pt->node_pt(9);
2253  }
2254 
2255  // Create new nodes
2256  // If we have an enriched element then don't need to construct centroid
2257  // node
2258  if (nnode == 15)
2259  {
2260  node0_pt = el_pt->node_pt(14);
2261  el1_pt->node_pt(2) = node0_pt;
2262  node5_pt = el1_pt->construct_node(11, time_stepper_pt);
2263  node6_pt = el1_pt->construct_node(13, time_stepper_pt);
2264  node9_pt = el1_pt->construct_node(12, time_stepper_pt);
2265 
2266  // Copy others over
2267  el1_pt->node_pt(10) = el_pt->node_pt(10);
2268  }
2269  // If not enriched we do
2270  else
2271  {
2272  node0_pt = el1_pt->construct_node(2, time_stepper_pt);
2273  }
2274  if (nnode_1d == 3)
2275  {
2276  node1_pt = el1_pt->construct_boundary_node(5, time_stepper_pt);
2277  node2_pt = el1_pt->construct_boundary_node(7, time_stepper_pt);
2278  node4_pt = el1_pt->construct_boundary_node(8, time_stepper_pt);
2279  }
2280 
2281 
2282  // Create second new element
2283  FiniteElement* el2_pt = new ELEMENT;
2284 
2285  // Copy existing nodes
2286  el2_pt->node_pt(0) = el_pt->node_pt(0);
2287  el2_pt->node_pt(1) = el_pt->node_pt(1);
2288  el2_pt->node_pt(2) = el_pt->node_pt(2);
2289  if (nnode_1d == 3)
2290  {
2291  el2_pt->node_pt(4) = el_pt->node_pt(4);
2292  el2_pt->node_pt(5) = el_pt->node_pt(5);
2293  el2_pt->node_pt(7) = el_pt->node_pt(7);
2294  }
2295 
2296  // Create new node
2297  if (nnode_1d == 3)
2298  {
2299  node3_pt = el2_pt->construct_boundary_node(8, time_stepper_pt);
2300  }
2301 
2302  // Copy existing new nodes
2303  el2_pt->node_pt(3) = node0_pt;
2304  if (nnode_1d == 3)
2305  {
2306  el2_pt->node_pt(6) = node1_pt;
2307  el2_pt->node_pt(9) = node2_pt;
2308  }
2309 
2310  // Copy and constuct other nodes for enriched elements
2311  if (nnode == 15)
2312  {
2313  el2_pt->node_pt(10) = node5_pt;
2314  el2_pt->node_pt(11) = el_pt->node_pt(11);
2315  node8_pt = el2_pt->construct_node(12, time_stepper_pt);
2316  node10_pt = el2_pt->construct_node(13, time_stepper_pt);
2317  }
2318 
2319  // Create third new element
2320  FiniteElement* el3_pt = new ELEMENT;
2321 
2322  // Copy existing nodes
2323  el3_pt->node_pt(1) = el_pt->node_pt(1);
2324  el3_pt->node_pt(2) = el_pt->node_pt(2);
2325  el3_pt->node_pt(3) = el_pt->node_pt(3);
2326  if (nnode_1d == 3)
2327  {
2328  el3_pt->node_pt(7) = el_pt->node_pt(7);
2329  el3_pt->node_pt(8) = el_pt->node_pt(8);
2330  el3_pt->node_pt(9) = el_pt->node_pt(9);
2331  }
2332 
2333  // Copy existing new nodes
2334  el3_pt->node_pt(0) = node0_pt;
2335  if (nnode_1d == 3)
2336  {
2337  el3_pt->node_pt(4) = node2_pt;
2338  el3_pt->node_pt(5) = node3_pt;
2339  el3_pt->node_pt(6) = node4_pt;
2340  }
2341 
2342  // Copy and constuct other nodes for enriched elements
2343  if (nnode == 15)
2344  {
2345  el3_pt->node_pt(10) = node6_pt;
2346  el3_pt->node_pt(11) = node10_pt;
2347  node7_pt = el3_pt->construct_node(12, time_stepper_pt);
2348  el3_pt->node_pt(13) = el_pt->node_pt(13);
2349  }
2350 
2351 
2352  // Create fourth new element
2353  FiniteElement* el4_pt = new ELEMENT;
2354 
2355  // Copy existing nodes
2356  el4_pt->node_pt(0) = el_pt->node_pt(0);
2357  el4_pt->node_pt(2) = el_pt->node_pt(2);
2358  el4_pt->node_pt(3) = el_pt->node_pt(3);
2359  if (nnode_1d == 3)
2360  {
2361  el4_pt->node_pt(5) = el_pt->node_pt(5);
2362  el4_pt->node_pt(6) = el_pt->node_pt(6);
2363  el4_pt->node_pt(8) = el_pt->node_pt(8);
2364  }
2365 
2366  // Copy existing new nodes
2367  el4_pt->node_pt(1) = node0_pt;
2368  if (nnode_1d == 3)
2369  {
2370  el4_pt->node_pt(4) = node1_pt;
2371  el4_pt->node_pt(7) = node3_pt;
2372  el4_pt->node_pt(9) = node4_pt;
2373  }
2374 
2375  // Copy all other nodes
2376  if (nnode == 15)
2377  {
2378  el4_pt->node_pt(10) = node9_pt;
2379  el4_pt->node_pt(11) = node8_pt;
2380  el4_pt->node_pt(12) = el_pt->node_pt(12);
2381  el4_pt->node_pt(13) = node7_pt;
2382  ;
2383  }
2384 
2385 
2386  // Add new elements and nodes
2387  new_or_retained_el_pt.push_back(el1_pt);
2388  new_or_retained_el_pt.push_back(el2_pt);
2389  new_or_retained_el_pt.push_back(el3_pt);
2390  new_or_retained_el_pt.push_back(el4_pt);
2391 
2392  // create a vector of the newly created elements
2393  Vector<FiniteElement*> temp_vec(4);
2394  temp_vec[0] = el1_pt;
2395  temp_vec[1] = el2_pt;
2396  temp_vec[2] = el3_pt;
2397  temp_vec[3] = el4_pt;
2398 
2399  // add the vector to the map
2400  old_to_new_corner_element_map.insert(
2401  std::pair<FiniteElement*, Vector<FiniteElement*>>(el_pt, temp_vec));
2402 
2403  if (nnode != 15)
2404  {
2405  this->add_node_pt(node0_pt);
2406  }
2407  this->add_node_pt(node1_pt);
2408  this->add_node_pt(node2_pt);
2409  this->add_node_pt(node3_pt);
2410  this->add_node_pt(node4_pt);
2411  if (nnode == 15)
2412  {
2413  this->add_node_pt(node5_pt);
2414  this->add_node_pt(node6_pt);
2415  this->add_node_pt(node7_pt);
2416  this->add_node_pt(node8_pt);
2417  this->add_node_pt(node9_pt);
2418  }
2419 
2420  // Set nodal positions
2421  for (unsigned i = 0; i < 3; i++)
2422  {
2423  // Only bother to set centroid if does not already exist
2424  if (nnode != 15)
2425  {
2426  node0_pt->x(i) =
2427  0.25 * (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i) +
2428  el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i));
2429  }
2430 
2431  if (nnode_1d == 3)
2432  {
2433  node1_pt->x(i) = 0.5 * (el_pt->node_pt(0)->x(i) + node0_pt->x(i));
2434  node2_pt->x(i) = 0.5 * (el_pt->node_pt(1)->x(i) + node0_pt->x(i));
2435  node3_pt->x(i) = 0.5 * (el_pt->node_pt(2)->x(i) + node0_pt->x(i));
2436  node4_pt->x(i) = 0.5 * (el_pt->node_pt(3)->x(i) + node0_pt->x(i));
2437  }
2438  }
2439 
2440 
2441  // Construct the four interior nodes if needed
2442  // and add to the list
2443  if (nnode == 15)
2444  {
2445  // Set the positions of the newly created mid-face nodes
2446  // New Node 5 lies in the plane between original nodes 0 1 centroid
2447  for (unsigned i = 0; i < 3; ++i)
2448  {
2449  node5_pt->x(i) =
2450  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i) +
2451  el_pt->node_pt(14)->x(i)) /
2452  3.0;
2453  }
2454 
2455  // New Node 6 lies in the plane between original nodes 1 3 centroid
2456  for (unsigned i = 0; i < 3; ++i)
2457  {
2458  node6_pt->x(i) =
2459  (el_pt->node_pt(1)->x(i) + el_pt->node_pt(3)->x(i) +
2460  el_pt->node_pt(14)->x(i)) /
2461  3.0;
2462  }
2463 
2464  // New Node 7 lies in the plane between original nodes 2 3 centroid
2465  for (unsigned i = 0; i < 3; ++i)
2466  {
2467  node7_pt->x(i) =
2468  (el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i) +
2469  el_pt->node_pt(14)->x(i)) /
2470  3.0;
2471  }
2472 
2473  // New Node 8 lies in the plane between original nodes 0 2 centroid
2474  for (unsigned i = 0; i < 3; ++i)
2475  {
2476  node8_pt->x(i) =
2477  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(2)->x(i) +
2478  el_pt->node_pt(14)->x(i)) /
2479  3.0;
2480  }
2481 
2482  // New Node 9 lies in the plane between original nodes 0 3 centroid
2483  for (unsigned i = 0; i < 3; ++i)
2484  {
2485  node9_pt->x(i) =
2486  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(3)->x(i) +
2487  el_pt->node_pt(14)->x(i)) /
2488  3.0;
2489  }
2490 
2491  // New Node 10 lies in the plane between original nodes 1 2 centroid
2492  for (unsigned i = 0; i < 3; ++i)
2493  {
2494  node10_pt->x(i) =
2495  (el_pt->node_pt(1)->x(i) + el_pt->node_pt(2)->x(i) +
2496  el_pt->node_pt(14)->x(i)) /
2497  3.0;
2498  }
2499 
2500  // Now create the new centroid nodes
2501 
2502  // First element
2503  Node* temp_nod_pt = el1_pt->construct_node(14, time_stepper_pt);
2504  for (unsigned i = 0; i < 3; ++i)
2505  {
2506  double av_pos = 0.0;
2507  for (unsigned j = 0; j < 4; j++)
2508  {
2509  av_pos += el1_pt->node_pt(j)->x(i);
2510  }
2511 
2512  temp_nod_pt->x(i) = 0.25 * av_pos;
2513  }
2514 
2515  this->add_node_pt(temp_nod_pt);
2516 
2517  // Second element
2518  temp_nod_pt = el2_pt->construct_node(14, time_stepper_pt);
2519  for (unsigned i = 0; i < 3; ++i)
2520  {
2521  double av_pos = 0.0;
2522  for (unsigned j = 0; j < 4; j++)
2523  {
2524  av_pos += el2_pt->node_pt(j)->x(i);
2525  }
2526  temp_nod_pt->x(i) = 0.25 * av_pos;
2527  }
2528  this->add_node_pt(temp_nod_pt);
2529 
2530  // Third element
2531  temp_nod_pt = el3_pt->construct_node(14, time_stepper_pt);
2532  for (unsigned i = 0; i < 3; ++i)
2533  {
2534  double av_pos = 0.0;
2535  for (unsigned j = 0; j < 4; j++)
2536  {
2537  av_pos += el3_pt->node_pt(j)->x(i);
2538  }
2539  temp_nod_pt->x(i) = 0.25 * av_pos;
2540  }
2541  this->add_node_pt(temp_nod_pt);
2542 
2543  // Fourth element
2544  temp_nod_pt = el4_pt->construct_node(14, time_stepper_pt);
2545  for (unsigned i = 0; i < 3; ++i)
2546  {
2547  double av_pos = 0.0;
2548  for (unsigned j = 0; j < 4; j++)
2549  {
2550  av_pos += el4_pt->node_pt(j)->x(i);
2551  }
2552  temp_nod_pt->x(i) = 0.25 * av_pos;
2553  }
2554  this->add_node_pt(temp_nod_pt);
2555  }
2556 
2557  // Kill the old element
2558  delete el_pt;
2559  }
2560  }
2561 
2562  // Flush element storage
2563  Element_pt.clear();
2564 
2565  // Copy across
2566  nel = new_or_retained_el_pt.size();
2567  Element_pt.resize(nel);
2568  for (unsigned e = 0; e < nel; e++)
2569  {
2570  Element_pt[e] = new_or_retained_el_pt[e];
2571  }
2572 
2573  // Setup boundary lookup scheme again
2575 
2576  // -------------------------------------------------------------------------
2577  // The various boundary/region lookups now need updating to account for the
2578  // newly added/removed elements. This will be done in two stages:
2579  // Step 1: Add the new elements to the vector of elements in the same region
2580  // as the original corner element, and then delete the originals.
2581  // Updating this lookup makes things easier in the following step.
2582  // Step 2: Regenerate the two more specific lookups: One which gives the
2583  // elements on a given boundary in a given region, and the other
2584  // which maps elements on a given boundary in a given region to the
2585  // element's face index on that boundary.
2586  //
2587  // N.B. the lookup Triangular_facet_vertex_boundary_coordinate is setup in
2588  // the call to setup_boundary_element_info() above so doesn't need
2589  // additional work.
2590 
2591  // if we have no regions then we have no lookups to update so we're done
2592  // here
2593  if (Region_attribute.size() == 0)
2594  {
2595  return;
2596  }
2597  // if we haven't had to split any corner elements then don't need to fiddle
2598  // with the lookups
2599  if (old_to_new_corner_element_map.size() == 0)
2600  {
2601  oomph_info << "\nNo corner elements need splitting\n\n";
2602  return;
2603  }
2604 
2605  // ------------------------------------------
2606  // Step 1: update the region element lookup
2607 
2608  // loop over the map of old corner elements which have been split
2609  for (std::map<FiniteElement*, Vector<FiniteElement*>>::iterator map_it =
2610  old_to_new_corner_element_map.begin();
2611  map_it != old_to_new_corner_element_map.end();
2612  map_it++)
2613  {
2614  // extract the old and new elements from the map
2615  FiniteElement* original_el_pt = map_it->first;
2616  Vector<FiniteElement*> new_el_pt = map_it->second;
2617 
2618  unsigned original_region_index = 0;
2619 
2620 #ifdef PARANOID
2621  // flag for paranoia, if for some reason we don't find the original corner
2622  // element in any of the regions
2623  bool found = false;
2624 #endif
2625 
2626  Vector<FiniteElement*>::iterator region_element_it;
2627 
2628  // loop over the regions and look for this original corner element to find
2629  // out which region it used to be in, so we can add the new elements to
2630  // the same region.
2631  for (unsigned region_index = 0; region_index < Region_element_pt.size();
2632  region_index++)
2633  {
2634  // for each region, search the vector of elements in this region for the
2635  // original corner element
2636  region_element_it = std::find(Region_element_pt[region_index].begin(),
2637  Region_element_pt[region_index].end(),
2638  original_el_pt);
2639 
2640  // if the iterator hasn't reached the end then we've found it
2641  if (region_element_it != Region_element_pt[region_index].end())
2642  {
2643  // save the region index we're at
2644  original_region_index = region_index;
2645 
2646 #ifdef PARANOID
2647  // update the paranoid flag
2648  found = true;
2649 #endif
2650 
2651  // regions can't overlap, so once we've found one we're done
2652  break;
2653  }
2654  }
2655 
2656 #ifdef PARANOID
2657  if (!found)
2658  {
2659  std::ostringstream error_message;
2660  error_message
2661  << "The corner element being split does not appear to be in any "
2662  << "region, so something has gone wrong with the region lookup "
2663  "scheme\n";
2664 
2665  throw OomphLibError(error_message.str(),
2666  OOMPH_CURRENT_FUNCTION,
2667  OOMPH_EXCEPTION_LOCATION);
2668  }
2669 #endif
2670 
2671  // Now update the basic region lookup. The iterator will still point to
2672  // the original corner element in the lookup, so we can delete this easily
2673  Region_element_pt[original_region_index].erase(region_element_it);
2674  for (unsigned i = 0; i < 4; i++)
2675  {
2676  Region_element_pt[original_region_index].push_back(new_el_pt[i]);
2677  }
2678  }
2679  // ------------------------------------------
2680  // Step 2: Clear and regenerate lookups
2681 
2684 
2687 
2688  for (unsigned b = 0; b < nboundary(); b++)
2689  {
2690  // Loop over elements next to that boundary
2691  unsigned nel = this->nboundary_element(b);
2692  for (unsigned e = 0; e < nel; e++)
2693  {
2694  FiniteElement* el_pt = boundary_element_pt(b, e);
2695 
2696  // now search for it in each region
2697  for (unsigned r_index = 0; r_index < Region_attribute.size(); r_index++)
2698  {
2699  unsigned region_id = static_cast<unsigned>(Region_attribute[r_index]);
2700 
2702  std::find(Region_element_pt[r_index].begin(),
2703  Region_element_pt[r_index].end(),
2704  el_pt);
2705 
2706  // if we find this element in the current region, then update our
2707  // lookups
2708  if (it != Region_element_pt[r_index].end())
2709  {
2710  Boundary_region_element_pt[b][region_id].push_back(el_pt);
2711 
2712  unsigned face_index = face_index_at_boundary(b, e);
2713  Face_index_region_at_boundary[b][region_id].push_back(face_index);
2714  }
2715  }
2716  }
2717  }
2718 
2719  oomph_info << "\nNumber of outer corner elements split: "
2720  << old_to_new_corner_element_map.size() << "\n\n";
2721 
2722  } // end split_elements_in_corners()
2723 } // namespace oomph
2724 
2725 #endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh....
Definition: nodes.h:2242
Base class for upgraded disk-like GeomObject (i.e. 2D surface in 3D space) with specification of boun...
Information for documentation of results: Directory and file number to enable output in the form RESL...
bool is_doc_enabled() const
Are we documenting?
void disable_doc()
Disable documentation.
std::string & label()
String used (e.g.) for labeling output files.
std::string directory() const
Output directory.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5018
A general Finite Element class.
Definition: elements.h:1317
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition: elements.h:2513
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2179
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2214
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2222
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition: elements.h:2542
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5123
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
A general mesh class.
Definition: mesh.h:67
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition: mesh.h:75
bool Lookup_for_elements_next_boundary_is_setup
Flag to indicate that the lookup schemes for elements that are adjacent to the boundaries has been se...
Definition: mesh.h:87
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:190
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
For the e-th finite element on boundary b, return int to indicate the face_index of the face adjacent...
Definition: mesh.h:896
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition: nodes.h:1373
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2394
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....
Definition: nodes.cc:2379
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: tet_mesh.h:847
virtual ~TetMeshBase()
Destructor (empty)
Definition: tet_mesh.h:859
void snap_nodes_onto_geometric_objects()
Move the nodes on boundaries with associated GeomObjects so that they exactly coincide with the geome...
Definition: tet_mesh.cc:483
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,...
Definition: tet_mesh.h:1645
int face_index_at_boundary_in_region(const unsigned &b, const unsigned &r, const unsigned &e) const
Return face index of the e-th element adjacent to boundary b in region r.
Definition: tet_mesh.h:1020
Vector< double > Region_attribute
Vector of attributes associated with the elements in each region NOTE: double is enforced on us by te...
Definition: tet_mesh.h:1213
double region_attribute(const unsigned &i)
Return the i-th region attribute (here only used as the (assumed to be unsigned) region id.
Definition: tet_mesh.h:1092
void setup_boundary_coordinates(const unsigned &b, const bool &switch_normal)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
Definition: tet_mesh.h:919
unsigned nregion()
Return the number of regions specified by attributes.
Definition: tet_mesh.h:1046
TimeStepper * Time_stepper_pt
Timestepper used to build nodes.
Definition: tet_mesh.h:1243
unsigned nboundary_element_in_region(const unsigned &b, const unsigned &r) const
Return the number of elements adjacent to boundary b in region r.
Definition: tet_mesh.h:983
std::map< unsigned, TetMeshFacetedSurface * > Tet_mesh_faceted_surface_pt
Reverse lookup scheme: Pointer to faceted surface (if any!) associated with boundary b.
Definition: tet_mesh.h:1231
Vector< TetMeshFacetedSurface * > Internal_surface_pt
Vector to faceted surfaces that define internal boundaries.
Definition: tet_mesh.h:1227
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
Storage for elements adjacent to a boundary in a particular region.
Definition: tet_mesh.h:1217
void operator=(const TetMeshBase &)=delete
Broken assignment operator.
std::map< unsigned, TetMeshFacet * > Tet_mesh_facet_pt
Reverse lookup scheme: Pointer to facet (if any!) associated with boundary b.
Definition: tet_mesh.h:1235
void setup_boundary_coordinates(const unsigned &b)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
Definition: tet_mesh.h:888
TetMeshFacetedClosedSurface * Outer_boundary_pt
Faceted surface that defines outer boundaries.
Definition: tet_mesh.h:1224
FiniteElement * boundary_element_in_region_pt(const unsigned &b, const unsigned &r, const unsigned &e) const
Return pointer to the e-th element adjacent to boundary b in region r.
Definition: tet_mesh.h:1002
static double Tolerance_for_boundary_finding
Global static data that specifies the permitted error in the setup of the boundary coordinates.
Definition: tet_mesh.h:863
void assess_mesh_quality(std::ofstream &some_file)
Assess mesh quality: Ratio of max. edge length to min. height, so if it's very large it's BAAAAAD.
Definition: tet_mesh.cc:373
FiniteElement * region_element_pt(const unsigned &r, const unsigned &e)
Return the e-th element in the r-th region.
Definition: tet_mesh.h:1098
TetMeshBase()
Constructor.
Definition: tet_mesh.h:850
void snap_to_quadratic_surface(const Vector< unsigned > &boundary_id, const std::string &quadratic_surface_file_name, const bool &switch_normal)
Snap boundaries specified by the IDs listed in boundary_id to a quadratric surface,...
Definition: tet_mesh.h:1165
void setup_boundary_element_info()
Setup lookup schemes which establish which elements are located next to mesh's boundaries (wrapper to...
Definition: tet_mesh.h:1193
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Storage for the face index adjacent to a boundary in a particular region.
Definition: tet_mesh.h:1221
std::map< unsigned, Vector< Vector< double > > > Triangular_facet_vertex_boundary_coordinate
Boundary coordinates of vertices in triangular facets associated with given boundary....
Definition: tet_mesh.h:1240
TetMeshBase(const TetMeshBase &node)=delete
Broken copy constructor.
Vector< Vector< FiniteElement * > > Region_element_pt
Vectors of vectors of elements in each region (note: this just stores them; the region IDs are contai...
Definition: tet_mesh.h:1208
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...
Definition: tet_mesh.h:2136
void setup_boundary_coordinates(const unsigned &b, std::ofstream &outfile)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
Definition: tet_mesh.h:974
unsigned nregion_element(const unsigned &r)
Return the number of elements in region r.
Definition: tet_mesh.h:1052
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:168
Vector< TetMeshVertex * > Vertex_pt
Pointer to vertices.
Definition: tet_mesh.h:280
void set_vertex_pt(const unsigned &j, TetMeshVertex *vertex_pt)
Set pointer to j-th vertex and pass one-based boundary id that may already have been set for this fac...
Definition: tet_mesh.h:194
TetMeshVertex * vertex_pt(const unsigned &j) const
Pointer to j-th vertex (const)
Definition: tet_mesh.h:187
void output(std::ostream &outfile) const
Output.
Definition: tet_mesh.h:263
void set_one_based_region_that_facet_is_embedded_in(const unsigned &one_based_region_id)
Facet is to be embedded in specified one-based region.
Definition: tet_mesh.h:249
std::set< unsigned > One_based_adjacent_region_id
Set of one-based adjacent region ids; no adjacent region if it's zero.
Definition: tet_mesh.h:287
bool facet_is_embedded_in_a_specified_region()
Boolean indicating that facet is embedded in a specified region.
Definition: tet_mesh.h:243
void set_one_based_boundary_id(const unsigned &one_based_id)
Set (one-based!) boundary IDs this facet lives on. Passed to any existing vertices and to any future ...
Definition: tet_mesh.h:214
unsigned one_based_boundary_id() const
Constant access to (one-based!) boundary IDs this facet lives on.
Definition: tet_mesh.h:207
unsigned nvertex() const
Number of vertices.
Definition: tet_mesh.h:181
unsigned One_based_boundary_id
(One-based!) boundary IDs this facet lives on
Definition: tet_mesh.h:283
unsigned One_based_region_id_that_facet_is_embedded_in
Facet is to be embedded in specified one-based region. Defaults to zero, indicating that its not embe...
Definition: tet_mesh.h:292
void set_one_based_adjacent_region_id(const unsigned &one_based_id)
Set (one-based!) region ID this facet is adjacent to. Specification of zero argument is harmless as i...
Definition: tet_mesh.h:231
unsigned one_based_region_that_facet_is_embedded_in()
Which (one-based) region is facet embedded in (if zero, it's not embedded in any region)
Definition: tet_mesh.h:257
std::set< unsigned > one_based_adjacent_region_id() const
Return set of (one-based!) region IDs this facet is adjacent to.
Definition: tet_mesh.h:237
TetMeshFacet(const unsigned &nvertex)
Constructor: Specify number of vertices.
Definition: tet_mesh.h:171
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: tet_mesh.h:820
virtual ~TetMeshFacetedClosedSurfaceForRemesh()
Destructor. Delete allocated memory.
Definition: tet_mesh.cc:73
TetMeshFacetedClosedSurfaceForRemesh(Vector< Node * > const &vertex_node_pt, Vector< Vector< unsigned >> const &facet_connectivity, Vector< unsigned > const &facet_boundary_id)
Constructor for a FacetedSurface created from a list of nodes and connectivity information....
Definition: tet_mesh.cc:40
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:724
TetMeshFacetedClosedSurface()
Constructor:
Definition: tet_mesh.h:727
bool Faceted_volume_represents_hole_for_gmsh
Does closed surface represent hole for gmsh?
Definition: tet_mesh.h:810
void enable_faceted_volume_represents_hole_for_gmsh()
Declare closed surface to represent hole for gmsh.
Definition: tet_mesh.h:736
const int & region_id_for_tetgen(const unsigned &j) const
Return the (zero-based) region ID of j-th internal point for tetgen. Negative if it's actually a hole...
Definition: tet_mesh.h:784
bool faceted_volume_represents_hole_for_gmsh() const
Does closed surface represent hole for gmsh?
Definition: tet_mesh.h:748
Vector< std::pair< Vector< double >, int > > Internal_point_for_tetgen
Storage for internal points for tetgen. Stores pair of: – Vector containing coordinates of internal p...
Definition: tet_mesh.h:807
virtual ~TetMeshFacetedClosedSurface()
Empty destructor.
Definition: tet_mesh.h:733
unsigned ninternal_point_for_tetgen()
Number of internal points (identifying either regions or holes) for tetgen.
Definition: tet_mesh.h:777
bool internal_point_identifies_region_for_tetgen(const unsigned &j)
Is j-th internal point for tetgen associated with a region?
Definition: tet_mesh.h:797
bool internal_point_identifies_hole_for_tetgen(const unsigned &j)
Is j-th internal point for tetgen associated with a hole?
Definition: tet_mesh.h:791
const double & internal_point_for_tetgen(const unsigned &j, const unsigned &i) const
i=th coordinate of the j-th internal point for tetgen
Definition: tet_mesh.h:754
void disable_faceted_volume_represents_hole_for_gmsh()
Declare closed surface NOT to represent hole for gmsh.
Definition: tet_mesh.h:742
void set_region_for_tetgen(const unsigned &region_id, const Vector< double > &region_point)
Specify a region; pass (zero-based) region ID and coordinate of point in region for tetgen.
Definition: tet_mesh.h:768
void set_hole_for_tetgen(const Vector< double > &hole_point)
Specify coordinate of hole for tetgen.
Definition: tet_mesh.h:761
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:306
unsigned nvertex_on_facet(const unsigned &j) const
Number of vertices defining the j-th facet.
Definition: tet_mesh.h:349
bool Boundaries_can_be_split_in_tetgen
Boolean to indicate whether extra vertices can be added on the boundary in tetgen.
Definition: tet_mesh.h:676
void disable_boundaries_can_be_split_in_tetgen()
Test whether boundaries can be split in tetgen.
Definition: tet_mesh.h:367
virtual void boundary_zeta01(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:603
Vector< TetMeshVertex * > Vertex_pt
Vector pointers to vertices.
Definition: tet_mesh.h:669
unsigned nfacet() const
Number of facets.
Definition: tet_mesh.h:325
Vector< Vector< unsigned > > Facet_vertex_index_in_tetgen
Facet connectivity: Facet_vertex_index[f][j] is the index of the j-th vertex (in the Vertex_pt vector...
Definition: tet_mesh.h:680
void output_paraview(const std::string &filename) const
Outputs the faceted surface into a file with the specified name in the Paraview format....
Definition: tet_mesh.h:589
DiskLikeGeomObjectWithBoundaries * Geom_object_with_boundaries_pt
GeomObject with boundaries associated with this surface.
Definition: tet_mesh.h:683
void output_paraview(std::ostream &outfile) const
Outputs the faceted surface into a specified file in the Paraview format for viewing in Paraview....
Definition: tet_mesh.h:418
virtual void boundary_zeta12(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:622
bool boundaries_can_be_split_in_tetgen()
Test whether boundary can be split in tetgen.
Definition: tet_mesh.h:355
Vector< TetMeshFacet * > Facet_pt
Vector of pointers to facets.
Definition: tet_mesh.h:672
unsigned one_based_facet_boundary_id(const unsigned &j) const
One-based boundary id of j-th facet.
Definition: tet_mesh.h:331
TetMeshFacet * facet_pt(const unsigned &j) const
Pointer to j-th facet.
Definition: tet_mesh.h:373
Vector< unsigned > vertex_index_in_tetgen(const unsigned &f)
Facet connectivity: vertex_index[j] is the index of the j-th vertex (in the Vertex_pt vector) in face...
Definition: tet_mesh.h:658
void enable_boundaries_can_be_split_in_tetgen()
Test whether boundaries can be split in tetgen.
Definition: tet_mesh.h:361
virtual void boundary_zeta20(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:641
void setup_facet_connectivity_for_tetgen()
Setup facet connectivity for tetgen.
Definition: tet_mesh.h:688
double vertex_coordinate(const unsigned &j, const unsigned &i) const
i-th coordinate of j-th vertex
Definition: tet_mesh.h:343
unsigned nvertex() const
Number of vertices.
Definition: tet_mesh.h:319
TetMeshFacetedSurface()
Constructor:
Definition: tet_mesh.h:309
void output(const std::string &filename) const
Output.
Definition: tet_mesh.h:403
TetMeshVertex * vertex_pt(const unsigned &j) const
Pointer to j-th vertex.
Definition: tet_mesh.h:379
void output(std::ostream &outfile) const
Output.
Definition: tet_mesh.h:392
virtual ~TetMeshFacetedSurface()
Empty destructor.
Definition: tet_mesh.h:316
DiskLikeGeomObjectWithBoundaries * geom_object_with_boundaries_pt()
Access to GeomObject with boundaries associated with this surface (Null if there isn't one!...
Definition: tet_mesh.h:386
unsigned one_based_vertex_boundary_id(const unsigned &j) const
First (of possibly multiple) one-based boundary id of j-th vertex.
Definition: tet_mesh.h:337
Vertex for Tet mesh generation. Can lie on multiple boundaries (identified via one-based enumeration!...
Definition: tet_mesh.h:50
void set_zeta_in_geom_object(const Vector< double > &zeta)
Set intrinisic coordinates in GeomObject.
Definition: tet_mesh.h:97
Vector< double > X
Coordinate vector.
Definition: tet_mesh.h:147
TetMeshVertex(Node *const &node_pt)
Definition: tet_mesh.h:72
Vector< double > zeta_in_geom_object() const
Get intrinisic coordinates in GeomObject (returns zero sized vector if no such coordinates have been ...
Definition: tet_mesh.h:118
std::set< unsigned > One_based_boundary_id
Set of (one-based!) boundary IDs this vertex lives on.
Definition: tet_mesh.h:150
unsigned one_based_boundary_id() const
First (of possibly multiple) one-based boundary id.
Definition: tet_mesh.h:130
Vector< double > Zeta_in_geom_object
Intrinisic coordinates in GeomObject (zero sized if there isn't one.
Definition: tet_mesh.h:154
void set_one_based_boundary_id(const unsigned &id)
Set of (one-based!) boundary IDs this vertex lives on.
Definition: tet_mesh.h:141
double x(const unsigned &i) const
i-th coordinate
Definition: tet_mesh.h:124
TetMeshVertex(const Vector< double > &x)
Constructor: Pass coordinates (length 3!)
Definition: tet_mesh.h:56
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...