gmsh_tet_mesh.template.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-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 
27 #ifndef OOMPH_GMSH_TET_MESH_HEADER
28 #define OOMPH_GMSH_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 #include <algorithm>
36 #include <iterator>
37 
38 // OOMPH-LIB Headers
40 #include "../generic/sample_point_parameters.h"
41 #include "../generic/mesh_as_geometric_object.h"
42 #include "../generic/projection.h"
43 
44 namespace oomph
45 {
46  //=========================================================================
47  /// Class to collate parameters for Gmsh mesh generation
48  //=========================================================================
50  {
51  public:
52  /// Specify outer boundary of domain to be meshed.
53  /// Other parameters get default values and can be set via
54  /// member functions
55  GmshParameters(TetMeshFacetedClosedSurface* const& outer_boundary_pt,
56  const std::string& gmsh_command_line_invocation)
58  Element_volume(1.0),
59  Target_size_file_name(".gmsh_target_size_for_adaptation.dat"),
60  Geo_and_msh_file_stem(".geo_file"),
62  Max_element_size(1.0),
63  Min_element_size(0.01),
67  5),
73  {
74  }
75 
76  /// Outer boundary
77  TetMeshFacetedClosedSurface*& outer_boundary_pt()
78  {
79  return Outer_boundary_pt;
80  }
81 
82  /// Internal boundaries
83  Vector<TetMeshFacetedSurface*>& internal_surface_pt()
84  {
85  return Internal_surface_pt;
86  }
87 
88  /// Uniform target element volume
89  double& element_volume()
90  {
91  return Element_volume;
92  }
93 
94  /// Filename for target volumes (for system-call based transfer to
95  /// gmsh during mesh adaptation). Default:
96  /// .gmsh_target_size_for_adaptation.dat
97  std::string& target_size_file_name()
98  {
99  return Target_size_file_name;
100  }
101 
102  /// String to be issued via system command to activate gmsh
104  {
106  }
107 
108  /// Stem for geo and msh files (input/output to command-line gmsh
109  /// invocation)
110  std::string& geo_and_msh_file_stem()
111  {
112  return Geo_and_msh_file_stem;
113  }
114 
115 
116  /// Max. element size during refinement
118  {
119  return Max_element_size;
120  }
121 
122  /// Min. element size during refinement
124  {
125  return Min_element_size;
126  }
127 
128  /// Max. permitted edge ratio
130  {
132  }
133 
134  /// (Isotropic) grid spacing for target size transfer
136  {
138  }
139 
140  /// Target size is transferred onto regular grid (for gmsh) by
141  /// locate zeta. We try to find the exact point in the existing
142  /// mesh but if we fail to converge from the nearest specified number
143  /// of sample points we use the nearest of those.
145  {
147  }
148 
149  /// Stem for filename used to doc target element sizes on
150  /// gmsh grid. No doc if stem is equal to empty string (or counter
151  /// is negative)
153  {
155  }
156 
157  /// Counter for filename used to doc target element sizes on
158  /// gmsh grid. No doc if stem is equal to empty string (or counter
159  /// is negative)
161  {
163  }
164 
165  /// Is projection of old solution onto new mesh disabled?
167  {
168  return Projection_is_disabled;
169  }
170 
171  /// Disable projection of old solution onto new mesh
173  {
174  Projection_is_disabled = true;
175  }
176 
177  /// Disable projection of old solution onto new mesh
179  {
180  Projection_is_disabled = false;
181  }
182 
183  /// Output filename for gmsh on-screen output
185  {
187  }
188 
189  /// Counter for marker that indicates where we are
190  /// in gmsh on-screen output
192  {
194  }
195 
196  private:
197  /// Pointer to outer boundary
198  TetMeshFacetedClosedSurface* Outer_boundary_pt;
199 
200  /// Internal boundaries
201  Vector<TetMeshFacetedSurface*> Internal_surface_pt;
202 
203  /// Uniform element volume
205 
206  /// Filename for target volume (for system-call based transfer
207  /// to gmsh during mesh adaptation)
209 
210  /// Stem for geo and msh files (input/output to command-line gmsh
211  /// invocation)
213 
214  /// Gmsh command line invocation string
216 
217  /// Max. element size during refinement
219 
220  /// Min. element size during refinement
222 
223  /// Max edge ratio before remesh gets triggered
225 
226  /// (Isotropic) grid spacing for target size transfer
228 
229  /// Target size is transferred onto regular grid (for gmsh) by
230  /// locate zeta. We try to find the exact point in the existing
231  /// mesh but if we fail to converge from the nearest specified number of
232  /// sample points we use the nearest of those.
233  unsigned
235 
236  /// Stem for filename used to doc target element sizes on
237  /// gmsh grid. No doc if stem is equal to empty string (or counter
238  /// is negative)
240 
241  /// Counter for filename used to doc target element sizes on
242  /// gmsh grid. No doc if stem is equal to empty string (or counter
243  /// is negative)
245 
246  /// Is projection of old solution onto new mesh disabled?
248 
249  /// Output filename for gmsh on-screen output
251 
252  /// Counter for marker that indicates where we are
253  /// in gmsh on-screen output
255  };
256 
257  /// //////////////////////////////////////////////////////////////////
258  /// //////////////////////////////////////////////////////////////////
259  /// //////////////////////////////////////////////////////////////////
260 
261 
262  //=========================================================================
263  /// Helper class to keep track of edges in tet mesh generation
264  //=========================================================================
265  class TetEdge
266  {
267  public:
268  /// Constructor: Pass two vertices, identified by their indices
269  /// Edge "direction" is from lower vertex to higher vertex id so
270  /// can compare if we're dealing with the same one...
271  TetEdge(const unsigned& vertex1, const unsigned& vertex2)
272  {
273  if (vertex1 > vertex2)
274  {
275  Reversed = true;
276  Vertex_pair = std::make_pair(vertex2, vertex1);
277  }
278  else if (vertex1 < vertex2)
279  {
280  Reversed = false;
281  Vertex_pair = std::make_pair(vertex1, vertex2);
282  }
283  else
284  {
285  throw OomphLibError(
286  "Muppet! You can't build an edge from one vertex to the same vertex!",
287  OOMPH_CURRENT_FUNCTION,
288  OOMPH_EXCEPTION_LOCATION);
289  }
290  }
291 
292 
293  /// First vertex id
294  unsigned first_vertex_id() const
295  {
296  return Vertex_pair.first;
297  }
298 
299  /// Second vertex id
300  unsigned second_vertex_id() const
301  {
302  return Vertex_pair.second;
303  }
304 
305 
306  /// Edge is reversed in the sense that vertex1 actually has a higher
307  /// id than vertex2 (when specified in the constructor)
308  bool is_reversed() const
309  {
310  return Reversed;
311  }
312 
313  /// Comparison operator: Edges are identical if their sorted
314  /// (and therefore possibly reversed) vertex ids agree
315  bool operator==(const TetEdge& tet_edge) const
316  {
317  return ((tet_edge.first_vertex_id() == Vertex_pair.first) &&
318  (tet_edge.second_vertex_id() == Vertex_pair.second));
319  }
320 
321  /// Comparison operator. Lexicographic comparison based on
322  /// vertex ids
323  bool operator<(const TetEdge& tet_edge) const
324  {
325  if ((tet_edge.first_vertex_id() == Vertex_pair.first) &&
326  (tet_edge.second_vertex_id() == Vertex_pair.second))
327  {
328  return false;
329  }
330  else
331  {
332  if (tet_edge.first_vertex_id() == Vertex_pair.first)
333  {
334  return (tet_edge.second_vertex_id() < Vertex_pair.second);
335  }
336  else
337  {
338  return (tet_edge.first_vertex_id() < Vertex_pair.first);
339  }
340  }
341  }
342 
343 
344  private:
345  /// The vertices (sorted by vertex ids)
346  std::pair<unsigned, unsigned> Vertex_pair;
347 
348  /// Is it reversed? I.e. is the first input vertex stored after
349  /// the second one?
350  bool Reversed;
351  };
352 
353 
354  /// ///////////////////////////////////////////////////////////////////////
355  /// ///////////////////////////////////////////////////////////////////////
356  /// ///////////////////////////////////////////////////////////////////////
357 
358 
359  /// Forward declaration
360  template<class ELEMENT>
361  class GmshTetMesh;
362 
363 
364  //=========================================================================
365  // Gmsh-based Tet scaffold mesh
366  //=========================================================================
367  class GmshTetScaffoldMesh : public virtual TetMeshBase, public virtual Mesh
368  {
369  public:
370  /// We're friends with the actual mesh
371  template<class ELEMENT>
372  friend class GmshTetMesh;
373 
374  /// Build mesh, based on specified parameters. If boolean is set
375  /// to true, the target element sizes are read from file (used during
376  /// adaptation; otherwise uniform target size is used).
377  GmshTetScaffoldMesh(GmshParameters* gmsh_parameters_pt,
378  const bool& use_mesh_grading_from_file)
379  : Gmsh_parameters_pt(gmsh_parameters_pt)
380  {
381  double t_start = TimingHelpers::timer();
382 
383  // Create .geo file
384  write_geo_file(use_mesh_grading_from_file);
385 
386  oomph_info << "Time for writing geo file : "
387  << TimingHelpers::timer() - t_start << " sec " << std::endl;
388  t_start = TimingHelpers::timer();
389 
390  // Execute gmsh on command line
391  std::string gmsh_command_line_string = "";
392 
393  std::string gmsh_onscreen_output_file_name =
394  gmsh_parameters_pt->gmsh_onscreen_output_file_name();
395  std::ofstream gmsh_on_screen_output_file;
396  std::stringstream marker;
397  if (gmsh_onscreen_output_file_name != "")
398  {
399  marker << "\n\n====================================================\n"
400  << " gmsh invocation: "
401  << gmsh_parameters_pt->gmsh_onscreen_output_counter()
402  << std::endl
403  << "====================================================\n\n\n"
404  << std::endl;
405  gmsh_parameters_pt->gmsh_onscreen_output_counter()++;
406  oomph_info << marker.str();
407  gmsh_on_screen_output_file.open(gmsh_onscreen_output_file_name.c_str(),
408  std::ofstream::app);
409  gmsh_on_screen_output_file << marker.str();
410  gmsh_on_screen_output_file.close();
411  }
412 
413  gmsh_command_line_string +=
416  if (gmsh_onscreen_output_file_name != "")
417  {
418  gmsh_command_line_string += " >> " + gmsh_onscreen_output_file_name;
419  }
420 
421  // Note return flag isn't particularly well defined but we report it
422  // anyway to aid detection of problems...
423  int return_flag = system(gmsh_command_line_string.c_str());
424  oomph_info << "fyi: return from system command: " << return_flag
425  << std::endl;
426  oomph_info << "Time for gmsh system call : "
427  << TimingHelpers::timer() - t_start << " sec " << std::endl;
428  t_start = TimingHelpers::timer();
429 
430 
431  // Create the mesh
433 
434  oomph_info << "Time for creating mesh from msh file: "
435  << TimingHelpers::timer() - t_start << " sec " << std::endl;
436  }
437 
438  private:
439  /// Create mesh from msh file (created internally via disk-based operations)
441  {
442  // Create filename from stem
443  std::string mesh_file_name =
445 
446  // Keep around if we ever write a version where name is passed in.
447  /* // Check extension */
448  /* #ifdef PARANOID */
449  /* std::string mesh_file_stem= */
450  /* mesh_file_name.substr(0,mesh_file_name.length()-4); */
451  /* std::string test=mesh_file_stem+".msh"; */
452  /* if (test!=mesh_file_name) */
453  /* { */
454  /* std::string error_msg("msh file has wrong extension: "); */
455  /* error_msg += " " + mesh_file_name; */
456  /* throw OomphLibError(error_msg, */
457  /* OOMPH_CURRENT_FUNCTION, */
458  /* OOMPH_EXCEPTION_LOCATION); */
459  /* } */
460  /* #endif */
461 
462  // Open wide...
463  std::ifstream mesh_file(mesh_file_name.c_str(), std::ios_base::in);
464 
465 // Check that the file actually opened correctly
466 #ifdef PARANOID
467  if (!mesh_file.is_open())
468  {
469  std::string error_msg("Failed to open mesh file: ");
470  error_msg += "\"" + mesh_file_name + "\".";
471  throw OomphLibError(
472  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
473  }
474 #endif
475 
476  // First line: Must be "$MeshFormat"
477  //----------------------------------
478  std::string line;
479  mesh_file >> line;
480  if (line != "$MeshFormat")
481  {
482  std::string error_msg(
483  "First line has to contain the string \"$MeshFormat\"; ");
484  error_msg += " yours contains: " + line;
485  throw OomphLibError(
486  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
487  }
488 
489  // Rest of line
490  mesh_file >> line;
491 
492 
493  // Nodes
494  //------
495 
496  // Now keep reading until we find the nodes
497  bool keep_going = true;
498  while (keep_going)
499  {
500  std::string line;
501  mesh_file >> line;
502  if (line == "$Nodes")
503  {
504  keep_going = false;
505  }
506  }
507  unsigned nnod = 0;
508  mesh_file >> nnod;
509  std::map<unsigned, Vector<double>> node_coordinate;
510  for (unsigned j = 0; j < nnod; j++)
511  {
512  unsigned node_number = 0;
513  mesh_file >> node_number;
514 
515  // Read rest of line word by word
516  std::string s;
517  std::getline(mesh_file, s);
518  std::istringstream iss(s);
519  std::string sub;
520  while (iss >> sub)
521  {
522  node_coordinate[node_number].push_back(atof(sub.c_str()));
523  }
524  }
525 
526  mesh_file >> line;
527  if (line != "$EndNodes")
528  {
529  std::string error_msg("Line has to contain the string \"$EndNodes\"; ");
530  error_msg += " yours contains: " + line;
531  throw OomphLibError(
532  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
533  }
534 
535  // Rewind
536  mesh_file.clear();
537  mesh_file.seekg(0);
538 
539  mesh_file >> line;
540  if (line != "$MeshFormat")
541  {
542  std::string error_msg(
543  "First line has to contain the string \"$MeshFormat\"; ");
544  error_msg += " yours contains: " + line;
545  throw OomphLibError(
546  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
547  }
548 
549  // Rest of line
550  mesh_file >> line;
551 
552  // Storage for boundaries the nodes are on
553  std::map<unsigned, std::set<unsigned>> one_based_boundaries_of_node;
554 
555  // Elements
556  //---------
557 
558  // node number = boundary_node[one_based_bound_id][...]
559  std::map<unsigned, std::set<unsigned>> boundary_node;
560 
561  // element number = region_element[one_based_region_id][...]
562  std::map<unsigned, std::set<unsigned>> region_element;
563 
564  // Now keep reading until we find the nodes
565  keep_going = true;
566  while (keep_going)
567  {
568  std::string line;
569  mesh_file >> line;
570  if (line == "$Elements")
571  {
572  keep_going = false;
573  }
574  }
575 
576 
577  unsigned highest_one_based_boundary_id = 0;
578  unsigned n_tet_el = 0;
579 
580  unsigned nel = 0;
581  mesh_file >> nel;
582  std::map<unsigned, Vector<unsigned>> element_node_index;
583  for (unsigned e = 0; e < nel; e++)
584  {
585  // For each element the msh file provides:
586  //
587  // elm-number elm-type number-of-tags < tag > ... node-number-list
588  //
589  // By default, the first tag is the number of the physical entity to
590  // which the element belongs; the second is the number of the elementary
591  // geometrical entity to which the element belongs; the third is the
592  // number of mesh partitions to which the element belongs, followed by
593  // the partition ids (negative partition ids indicate ghost cells). A
594  // zero tag is equivalent to no tag. Gmsh and most codes using the MSH 2
595  // format require at least the first two tags (physical and elementary
596  // tags).
597  unsigned el_number = 0;
598  mesh_file >> el_number;
599 
600 
601  unsigned el_type = 0;
602  mesh_file >> el_type;
603 
604  switch (el_type)
605  {
606  case 1:
607  // oomph_info << "Line element" << std::endl;
608  break;
609 
610  case 2:
611  // oomph_info << "Triangle" << std::endl;
612  break;
613 
614  case 4:
615  n_tet_el++;
616  // oomph_info << "Tet" << std::endl;
617  break;
618 
619  case 15:
620  // oomph_info << "Point" << std::endl;
621  break;
622 
623  default:
624  std::string error_msg("Can't handle element type: ");
625  error_msg += oomph::StringConversion::to_string(el_type);
626  throw OomphLibError(
627  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
628  }
629 
630  unsigned ntags;
631  mesh_file >> ntags;
632 
633  // Gmesh produces at least two flags:
634 
635  // Physical tag = what we use as boundary or region IDs.
636  int physical_tag = 0;
637  mesh_file >> physical_tag;
638 
639  // Geometric tag; not something we use
640  int geom_tag = 0;
641  mesh_file >> geom_tag;
642 
643  Vector<int> other_tags;
644  for (unsigned i = 2; i < ntags; i++)
645  {
646  int tag = 0;
647  mesh_file >> tag;
648  other_tags.push_back(tag);
649  }
650 
651 
652  // Now read the rest: node numbers
653  // https://stackoverflow.com/questions/16991002/getting-a-line-from-a-file-and-then-reading-word-by-word-c
654  std::string s;
655  std::getline(mesh_file, s);
656  std::istringstream iss(s);
657  std::string sub;
658  Vector<int> other_ints;
659  while (iss >> sub)
660  {
661  other_ints.push_back(atoi(sub.c_str()));
662  }
663  unsigned n_el_nod = other_ints.size();
664  for (unsigned j = 0; j < n_el_nod; j++)
665  {
666  unsigned node_number = unsigned(other_ints[j]);
667  element_node_index[el_number].push_back(node_number);
668 
669  // If the element is a triangle, add node to boundary
670  // lookup scheme (if tag is not zero, i.e. hasn't been specified)
671  if (el_type == 2)
672  {
673  if (physical_tag != 0)
674  {
675  boundary_node[unsigned(physical_tag)].insert(node_number);
676  one_based_boundaries_of_node[node_number].insert(
677  unsigned(physical_tag));
678  if (unsigned(physical_tag) > highest_one_based_boundary_id)
679  {
680  highest_one_based_boundary_id = physical_tag;
681  }
682  }
683  }
684  }
685  // If it's a bulk element (tet) and the physical tag (region id)
686  // is not equal to 0 add it to region list
687  if (el_type == 4)
688  {
689  if (physical_tag != 0)
690  {
691  region_element[unsigned(physical_tag)].insert(n_tet_el);
692  }
693  }
694  }
695 
696  mesh_file >> line;
697  if (line != "$EndElements")
698  {
699  std::string error_msg(
700  "Line has to contain the string \"$EndElements\"; ");
701  error_msg += " yours contains: " + line;
702  throw OomphLibError(
703  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
704  }
705 
706  // Rewind
707  mesh_file.clear();
708  mesh_file.seekg(0);
709 
710  mesh_file >> line;
711  if (line != "$MeshFormat")
712  {
713  std::string error_msg(
714  "First line has to contain the string \"$MeshFormat\"; ");
715  error_msg += " yours contains: " + line;
716  throw OomphLibError(
717  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
718  }
719 
720  // Rest of line
721  mesh_file >> line;
722  mesh_file.close();
723 
724 
725  // Identify elements next to boundaries
726  //-------------------------------------
727 
728  // Now loop over tet elements and check if three their nodes are
729  // on given boundary
730  std::map<unsigned, std::set<unsigned>> element_next_to_boundary;
731  for (std::map<unsigned, Vector<unsigned>>::iterator it =
732  element_node_index.begin();
733  it != element_node_index.end();
734  it++)
735  {
736  unsigned el_number = (*it).first;
737  unsigned nnod = ((*it).second).size();
738  if (nnod == 4)
739  {
740  std::map<unsigned, unsigned> boundary_node_count;
741  for (unsigned j = 0; j < nnod; j++)
742  {
743  unsigned node_number = ((*it).second)[j];
744  std::map<unsigned, std::set<unsigned>>::iterator itt =
745  one_based_boundaries_of_node.find(node_number);
746  if (itt != one_based_boundaries_of_node.end())
747  {
748  for (std::set<unsigned>::iterator ittt = (itt->second).begin();
749  ittt != (itt->second).end();
750  ittt++)
751  {
752  unsigned one_based_boundary_id = (*ittt);
753  boundary_node_count[one_based_boundary_id]++;
754  }
755  }
756  }
757  for (std::map<unsigned, unsigned>::iterator itt =
758  boundary_node_count.begin();
759  itt != boundary_node_count.end();
760  itt++)
761  {
762  if ((*itt).second == 3)
763  {
764  element_next_to_boundary[(*itt).first].insert(el_number);
765  }
766  }
767  }
768  }
769 
770  this->set_nboundary(highest_one_based_boundary_id);
771 
772  // Done reading/processing; now move across
773  // ----------------------------------------
774 
775  // Translate enumeration
776  Vector<unsigned> oomph_lib_node_number(4);
777  oomph_lib_node_number[0] = 0;
778  oomph_lib_node_number[1] = 3;
779  oomph_lib_node_number[2] = 2;
780  oomph_lib_node_number[3] = 1;
781 
782  // make space to accomodate all this in the actual mesh
783  Element_pt.resize(n_tet_el);
784  Node_pt.resize(nnod);
785 
786  // Existing nodes
787  Vector<Node*> existing_node_pt(nnod, 0);
788 
789  // Now process the simplex tets only (they have four nodes!)
790  unsigned e = 0;
791  for (std::map<unsigned, Vector<unsigned>>::iterator it =
792  element_node_index.begin();
793  it != element_node_index.end();
794  it++)
795  {
796  unsigned el_nnod = (*it).second.size();
797  if (el_nnod == 4)
798  {
799  // Here comes the new element
800  TElement<3, 2>* el_pt = new TElement<3, 2>;
801 
802  // Store it
803  Element_pt[e] = el_pt;
804  e++;
805 
806  // Make/get nodes
807  for (unsigned j = 0; j < el_nnod; j++)
808  {
809  // Get global, one-based node number
810  unsigned node_number = (*it).second[j];
811 
812  // Does it already exist?
813  if (existing_node_pt[node_number - 1] != 0)
814  {
815  el_pt->node_pt(oomph_lib_node_number[j]) =
816  existing_node_pt[node_number - 1];
817  }
818  // Make new node
819  else
820  {
821  Node* nod_pt = 0;
822 
823  // Is it on the boundary?
824  if (one_based_boundaries_of_node[node_number].size() == 0)
825  {
826  // Make normal node
827  nod_pt = el_pt->construct_node(oomph_lib_node_number[j]);
828  Node_pt[node_number - 1] = nod_pt;
829  existing_node_pt[node_number - 1] = nod_pt;
830  }
831  // Make boundary node
832  else
833  {
834  nod_pt =
835  el_pt->construct_boundary_node(oomph_lib_node_number[j]);
836  Node_pt[node_number - 1] = nod_pt;
837  existing_node_pt[node_number - 1] = nod_pt;
838 
839  // Add to boundary lookup scheme
840  for (std::set<unsigned>::iterator it =
841  one_based_boundaries_of_node[node_number].begin();
842  it != one_based_boundaries_of_node[node_number].end();
843  it++)
844  {
845  add_boundary_node((*it) - 1, nod_pt);
846  }
847  }
848 
849  // Assign coordinates of new node
850  Vector<double> x(node_coordinate[node_number]);
851  for (unsigned i = 0; i < 3; i++)
852  {
853  nod_pt->x(i) = x[i];
854  }
855  }
856  }
857  } // end for tet element
858 
859  } // End of loop over all elements
860 
861 
862  // Setup region info. This is ugly because we're using
863  // a lookup scheme that was originally designed for tetgen...
864  unsigned n_region = region_element.size();
865  this->Region_element_pt.resize(n_region);
866  this->Region_attribute.resize(n_region);
867  unsigned region_count = 0;
868  for (std::map<unsigned, std::set<unsigned>>::iterator it =
869  region_element.begin();
870  it != region_element.end();
871  it++)
872  {
873  this->Region_element_pt[region_count].resize((*it).second.size());
874  unsigned one_based_region_id = (*it).first;
875  this->Region_attribute[region_count] = one_based_region_id - 1;
876  unsigned count = 0;
877  for (std::set<unsigned>::iterator itt = (*it).second.begin();
878  itt != (*it).second.end();
879  itt++)
880  {
881  unsigned one_based_el_number = (*itt);
882  this->Region_element_pt[region_count][count] =
883  dynamic_cast<FiniteElement*>(Element_pt[one_based_el_number - 1]);
884  count++;
885  }
886  region_count++;
887  }
888 
889 
890  // Keep alive for debugging
891  bool plot_for_debugging = false;
892  if (plot_for_debugging)
893  {
894  std::ofstream outfile;
895  std::string outfile_name;
896  std::string mesh_file_stem = "shite";
897 
898  // Output element nodes
899  //----------------------
900  outfile_name = mesh_file_stem + "_element_nodes.dat";
901  outfile.open(outfile_name.c_str());
902  for (std::map<unsigned, Vector<unsigned>>::iterator it =
903  element_node_index.begin();
904  it != element_node_index.end();
905  it++)
906  {
907  std::set<unsigned> shite_nodes;
908  unsigned el_id = (*it).first;
909  std::stringstream tmp_out;
910  for (Vector<unsigned>::iterator itt = (*it).second.begin();
911  itt != (*it).second.end();
912  itt++)
913  {
914  unsigned node_number = (*itt);
915  shite_nodes.insert(node_number);
916  Vector<double> x(node_coordinate[node_number]);
917  unsigned n = x.size();
918  for (unsigned i = 0; i < n; i++)
919  {
920  tmp_out << x[i] << " ";
921  }
922  tmp_out << node_number << std::endl;
923  }
924  std::string prefix = ", N=4, E=1, F=FEPOINT, ET=TETRAHEDRON";
925  std::string postfix = "1 2 3 4";
926  if ((*it).second.size() == 3)
927  {
928  prefix = ", N=3, E=1, F=FEPOINT, ET=TRIANGLE";
929  postfix = "1 2 3";
930  }
931  outfile << "ZONE T=\"one-based element id=" << el_id << "\"" << prefix
932  << std::endl;
933  outfile << tmp_out.str();
934  outfile << postfix << std::endl;
935  }
936  outfile.close();
937 
938  // Output boundary nodes
939  //----------------------
940  outfile_name = mesh_file_stem + "_boundary_nodes.dat";
941  outfile.open(outfile_name.c_str());
942  for (std::map<unsigned, std::set<unsigned>>::iterator it =
943  boundary_node.begin();
944  it != boundary_node.end();
945  it++)
946  {
947  unsigned one_based_boundary_id = (*it).first;
948  outfile << "ZONE T=\"one-based boundary id=" << one_based_boundary_id
949  << "\"" << std::endl;
950  for (std::set<unsigned>::iterator itt = (*it).second.begin();
951  itt != (*it).second.end();
952  itt++)
953  {
954  unsigned node_number = (*itt);
955  Vector<double> x(node_coordinate[node_number]);
956  unsigned n = x.size();
957  for (unsigned i = 0; i < n; i++)
958  {
959  outfile << x[i] << " ";
960  }
961  outfile << node_number << std::endl;
962  }
963  }
964  outfile.close();
965 
966 
967  // Output elements next to boundaries
968  //-----------------------------------
969  for (std::map<unsigned, std::set<unsigned>>::iterator it =
970  element_next_to_boundary.begin();
971  it != element_next_to_boundary.end();
972  it++)
973  {
974  unsigned one_based_boundary_id = (*it).first;
975  outfile_name =
976  mesh_file_stem + "_elements_next_to_boundary_" +
977  oomph::StringConversion::to_string(one_based_boundary_id) + ".dat";
978  outfile.open(outfile_name.c_str());
979  for (std::set<unsigned>::iterator itt = (*it).second.begin();
980  itt != (*it).second.end();
981  itt++)
982  {
983  outfile << "ZONE T=\"one-based boundary " << one_based_boundary_id
984  << "\", N=4, E=1, F=FEPOINT, ET=TETRAHEDRON" << std::endl;
985  unsigned el_number = (*itt);
986  unsigned nnod = element_node_index[el_number].size();
987  for (unsigned j = 0; j < nnod; j++)
988  {
989  unsigned node_number = element_node_index[el_number][j];
990  Vector<double> x(node_coordinate[node_number]);
991  unsigned n = x.size();
992  for (unsigned i = 0; i < n; i++)
993  {
994  outfile << x[i] << " ";
995  }
996  outfile << std::endl;
997  }
998  outfile << "1 2 3 4" << std::endl;
999  }
1000  outfile.close();
1001  }
1002 
1003 
1004  // Compute/report total volume for sanity check
1005  {
1006  double vol = 0.0;
1007  unsigned nel = this->nelement();
1008  for (unsigned e = 0; e < nel; e++)
1009  {
1010  vol += this->finite_element_pt(e)->size();
1011  }
1012  oomph_info << "Total volume of all elements in scaffold mesh: " << vol
1013  << std::endl;
1014  }
1015  }
1016  }
1017 
1018 
1019  /// Write geo file for gmsh
1020  void write_geo_file(const bool& use_mesh_grading_from_file)
1021  {
1022  // Transfer data from parameters
1023  TetMeshFacetedClosedSurface* outer_boundary_pt =
1025 
1026  Vector<TetMeshFacetedSurface*> internal_surface_pt =
1028 
1029  double element_volume = Gmsh_parameters_pt->element_volume();
1030 
1031  std::string& target_size_file_name =
1033 
1034  // Write gmsh geo file:
1035  std::string filename =
1037  std::ofstream geo_file;
1038  geo_file.open(filename.c_str());
1039 
1040  geo_file << "// Uniform element size" << std::endl;
1041  geo_file << "//---------------------" << std::endl;
1042  geo_file << "lc=" << pow(element_volume, 1.0 / 3.0) << ";" << std::endl;
1043  geo_file << std::endl;
1044 
1045  // Outer boundary
1046  //===============
1047 
1048  // Create vertices
1049  geo_file << "// Outer box" << std::endl;
1050  geo_file << "//==========" << std::endl;
1051  geo_file << std::endl;
1052  geo_file << "// Vertices" << std::endl;
1053  geo_file << "//---------" << std::endl;
1054  std::map<TetMeshVertex*, unsigned> vertex_number;
1055  unsigned nv = outer_boundary_pt->nvertex();
1056  for (unsigned j = 0; j < nv; j++)
1057  {
1058  TetMeshVertex* vertex_pt = outer_boundary_pt->vertex_pt(j);
1059  vertex_number[vertex_pt] = j;
1060  geo_file << "Point(" << j + 1 << ")={";
1061  for (unsigned i = 0; i < 3; i++)
1062  {
1063  geo_file << vertex_pt->x(i) << ",";
1064  }
1065  geo_file << "lc};" << std::endl;
1066  }
1067 
1068 
1069  // Determine unique edges
1070  std::set<TetEdge> tet_edge_set;
1071  unsigned nfacet = outer_boundary_pt->nfacet();
1072  for (unsigned f = 0; f < nfacet; f++)
1073  {
1074  TetMeshFacet* facet_pt = outer_boundary_pt->facet_pt(f);
1075  unsigned nv = facet_pt->nvertex();
1076  for (unsigned j = 0; j < nv - 1; j++)
1077  {
1078  TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(j);
1079  TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(j + 1);
1080  TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1081  vertex_number[second_vertex_pt] + 1);
1082  tet_edge_set.insert(my_tet_edge);
1083  }
1084  TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(nv - 1);
1085  TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(0);
1086  TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1087  vertex_number[second_vertex_pt] + 1);
1088  tet_edge_set.insert(my_tet_edge);
1089  }
1090 
1091 
1092  geo_file << std::endl;
1093  geo_file << "// Edge of outer box" << std::endl;
1094  geo_file << "//------------------" << std::endl;
1095  unsigned count = 0;
1096  std::map<TetEdge, unsigned> tet_edge;
1097  for (std::set<TetEdge>::iterator it = tet_edge_set.begin();
1098  it != tet_edge_set.end();
1099  it++)
1100  {
1101  tet_edge.insert(std::make_pair((*it), count));
1102  geo_file << "Line(" << count + 1 << ")={" << (*it).first_vertex_id()
1103  << "," << (*it).second_vertex_id() << "};" << std::endl;
1104 
1105  count++;
1106  }
1107 
1108 
1109  geo_file << std::endl;
1110  geo_file << "// Faces of outer box" << std::endl;
1111  geo_file << "//-------------------" << std::endl;
1112  for (unsigned f = 0; f < nfacet; f++)
1113  {
1114  geo_file << "Line Loop(" << f + 1 << ")={";
1115 
1116  TetMeshFacet* facet_pt = outer_boundary_pt->facet_pt(f);
1117  unsigned nv = facet_pt->nvertex();
1118  for (unsigned j = 0; j < nv - 1; j++)
1119  {
1120  TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(j);
1121  TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(j + 1);
1122  TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1123  vertex_number[second_vertex_pt] + 1);
1124 
1125  std::map<TetEdge, unsigned>::iterator it = tet_edge.find(my_tet_edge);
1126  if (my_tet_edge.is_reversed())
1127  {
1128  geo_file << -int((it->second) + 1) << ",";
1129  }
1130  else
1131  {
1132  geo_file << ((it->second) + 1) << ",";
1133  }
1134  }
1135 
1136  TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(nv - 1);
1137  TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(0);
1138  TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1139  vertex_number[second_vertex_pt] + 1);
1140 
1141  std::map<TetEdge, unsigned>::iterator it = tet_edge.find(my_tet_edge);
1142  if (my_tet_edge.is_reversed())
1143  {
1144  geo_file << -int((it->second) + 1) << "};" << std::endl;
1145  }
1146  else
1147  {
1148  geo_file << ((it->second) + 1) << "};" << std::endl;
1149  }
1150  geo_file << "Plane Surface(" << f + 1 << ")={" << f + 1 << "};"
1151  << std::endl;
1152  }
1153 
1154 
1155  geo_file << std::endl;
1156  geo_file << "// Define Plane Surfaces bounding the volume" << std::endl;
1157  geo_file << "//------------------------------------------" << std::endl;
1158  geo_file << "Surface Loop(1) = {";
1159  for (unsigned f = 0; f < nfacet - 1; f++)
1160  {
1161  geo_file << f + 1 << ",";
1162  }
1163  geo_file << nfacet << "};" << std::endl;
1164 
1165 
1166  geo_file << std::endl;
1167  geo_file << "// Define one-based boundary IDs" << std::endl;
1168  geo_file << "//------------------------------" << std::endl;
1169  for (unsigned f = 0; f < nfacet; f++)
1170  {
1171  unsigned one_based_boundary_id =
1172  outer_boundary_pt->one_based_facet_boundary_id(f);
1173  geo_file << "Physical Surface(" << one_based_boundary_id << ") = {"
1174  << f + 1 << "};" << std::endl;
1175  }
1176 
1177 
1178  // Offsets before we start adding the various internal volumes/surfaces
1179  Vector<unsigned> volume_id_to_be_subtracted_off;
1180  unsigned nvertex_offset = outer_boundary_pt->nvertex();
1181  unsigned nfacet_offset = outer_boundary_pt->nfacet();
1182  unsigned nvolume_offset = 1;
1183  unsigned nline_offset = tet_edge.size();
1184 
1185 
1186  // Storage for one-based ids of surfaces to be embedded in
1187  // main volume
1188  Vector<unsigned> surfaces_to_be_embedded_in_main_volume;
1189  std::map<unsigned, Vector<unsigned>>
1190  surfaces_to_be_embedded_in_specified_one_based_region;
1191 
1192  // Now deal with internal faceted objects
1193  //=======================================
1194  unsigned n_internal = internal_surface_pt.size();
1195  for (unsigned i_internal = 0; i_internal < n_internal; i_internal++)
1196  {
1197  // Closed surface?
1198  TetMeshFacetedClosedSurface* closed_srf_pt =
1199  dynamic_cast<TetMeshFacetedClosedSurface*>(
1200  internal_surface_pt[i_internal]);
1201  bool inner_surface_is_closed = true;
1202  if (closed_srf_pt == 0)
1203  {
1204  inner_surface_is_closed = false;
1205  }
1206 
1207  // What it says
1208  unsigned number_of_volumes_created_for_this_internal_object = 0;
1209 
1210  // Create vertices
1211  geo_file << std::endl;
1212  geo_file << std::endl;
1213  geo_file << "// Inner faceted surface " << i_internal << std::endl;
1214  geo_file << "//==========================" << std::endl;
1215  geo_file << std::endl;
1216  geo_file << "// Vertices" << std::endl;
1217  geo_file << "//---------" << std::endl;
1218  std::map<TetMeshVertex*, unsigned> vertex_number;
1219  unsigned nv = internal_surface_pt[i_internal]->nvertex();
1220  for (unsigned j = 0; j < nv; j++)
1221  {
1222  TetMeshVertex* vertex_pt =
1223  internal_surface_pt[i_internal]->vertex_pt(j);
1224  vertex_number[vertex_pt] = nvertex_offset + j;
1225  geo_file << "Point(" << nvertex_offset + j + 1 << ")={";
1226  for (unsigned i = 0; i < 3; i++)
1227  {
1228  geo_file << vertex_pt->x(i) << ",";
1229  }
1230  geo_file << "lc};" << std::endl;
1231  }
1232 
1233  // Determine unique edges
1234  std::set<TetEdge> tet_edge_set;
1235  unsigned nfacet = internal_surface_pt[i_internal]->nfacet();
1236  for (unsigned f = 0; f < nfacet; f++)
1237  {
1238  TetMeshFacet* facet_pt = internal_surface_pt[i_internal]->facet_pt(f);
1239  unsigned nv = facet_pt->nvertex();
1240  for (unsigned j = 0; j < nv - 1; j++)
1241  {
1242  TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(j);
1243  TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(j + 1);
1244  TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1245  vertex_number[second_vertex_pt] + 1);
1246  tet_edge_set.insert(my_tet_edge);
1247  }
1248  TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(nv - 1);
1249  TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(0);
1250  TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1251  vertex_number[second_vertex_pt] + 1);
1252  tet_edge_set.insert(my_tet_edge);
1253  }
1254 
1255  geo_file << std::endl;
1256  geo_file << "// Edge of inner faceted surface" << std::endl;
1257  geo_file << "//------------------------------" << std::endl;
1258  unsigned count = 0;
1259  std::map<TetEdge, unsigned> tet_edge;
1260  for (std::set<TetEdge>::iterator it = tet_edge_set.begin();
1261  it != tet_edge_set.end();
1262  it++)
1263  {
1264  tet_edge.insert(std::make_pair((*it), nline_offset + count));
1265  geo_file << "Line(" << nline_offset + count + 1 << ")={"
1266  << (*it).first_vertex_id() << "," << (*it).second_vertex_id()
1267  << "};" << std::endl;
1268 
1269  count++;
1270  }
1271 
1272 
1273  geo_file << std::endl;
1274  geo_file << "// Faces of inner faceted surface" << std::endl;
1275  geo_file << "//-------------------------------" << std::endl;
1276  for (unsigned f = 0; f < nfacet; f++)
1277  {
1278  geo_file << "Line Loop(" << nfacet_offset + f + 1 << ")={";
1279 
1280  TetMeshFacet* facet_pt = internal_surface_pt[i_internal]->facet_pt(f);
1281  unsigned nv = facet_pt->nvertex();
1282  for (unsigned j = 0; j < nv - 1; j++)
1283  {
1284  TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(j);
1285  TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(j + 1);
1286  TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1287  vertex_number[second_vertex_pt] + 1);
1288 
1289  std::map<TetEdge, unsigned>::iterator it =
1290  tet_edge.find(my_tet_edge);
1291  if (my_tet_edge.is_reversed())
1292  {
1293  geo_file << -int((it->second) + 1) << ",";
1294  }
1295  else
1296  {
1297  geo_file << ((it->second) + 1) << ",";
1298  }
1299  }
1300 
1301  TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(nv - 1);
1302  TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(0);
1303  TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1304  vertex_number[second_vertex_pt] + 1);
1305 
1306  std::map<TetEdge, unsigned>::iterator it = tet_edge.find(my_tet_edge);
1307  if (my_tet_edge.is_reversed())
1308  {
1309  geo_file << -int((it->second) + 1) << "};" << std::endl;
1310  }
1311  else
1312  {
1313  geo_file << ((it->second) + 1) << "};" << std::endl;
1314  }
1315  geo_file << "Plane Surface(" << nfacet_offset + f + 1 << ")={"
1316  << nfacet_offset + f + 1 << "};" << std::endl;
1317 
1318  // Keep track of plane surfaces that we've created so we
1319  // can embed them in volume for non-closed surfaces
1320  bool facet_is_embedded_in_a_volume =
1321  facet_pt->facet_is_embedded_in_a_specified_region();
1322  if (facet_is_embedded_in_a_volume)
1323  {
1324  unsigned one_based_region_id =
1325  facet_pt->one_based_region_that_facet_is_embedded_in();
1326  if (one_based_region_id == 1)
1327  {
1328  surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset +
1329  f + 1);
1330  }
1331  else
1332  {
1333  surfaces_to_be_embedded_in_specified_one_based_region
1334  [one_based_region_id]
1335  .push_back(nfacet_offset + f + 1);
1336  }
1337  }
1338  else
1339  {
1340  // By default put all surfaces into main volume
1341  if (!inner_surface_is_closed)
1342  {
1343  surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset +
1344  f + 1);
1345  }
1346  }
1347  }
1348 
1349  // Store all region IDs defined by bounding facets
1350  std::set<unsigned> all_regions_id;
1351 
1352  // region_bounding_facet[r][i] returns the facet id (in gmsh
1353  // counting) of i-th facet bounding region r
1354  std::map<unsigned, Vector<unsigned>> region_bounding_facet;
1355 
1356  // outer_bounding_facet[i] returns the facet id (in gmsh
1357  // counting) that encloses the actual regions and acts as a
1358  // hole in the main volume (volume 1)
1359  Vector<unsigned> outer_bounding_facet;
1360 
1361  // Loop pver all facets to figure out which regions are bounded
1362  // by them
1363  for (unsigned f = 0; f < nfacet; f++)
1364  {
1365  TetMeshFacet* facet_pt = internal_surface_pt[i_internal]->facet_pt(f);
1366  std::set<unsigned> region_id(
1367  facet_pt->one_based_adjacent_region_id());
1368  unsigned nr = region_id.size();
1369  if (nr == 1) outer_bounding_facet.push_back(nfacet_offset + f + 1);
1370  if ((nr == 0) && inner_surface_is_closed)
1371  {
1372  // Add to list of plane surfaces that don't bound regions so we
1373  // can embed them in volume for non-closed surfaces
1374  surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset + f +
1375  1);
1376  }
1377  for (std::set<unsigned>::iterator it = region_id.begin();
1378  it != region_id.end();
1379  it++)
1380  {
1381  all_regions_id.insert((*it));
1382  region_bounding_facet[(*it)].push_back(nfacet_offset + f + 1);
1383  }
1384  }
1385 
1386  // Number of regions bounded by facets
1387  unsigned n_regions_bounded_by_facets = all_regions_id.size();
1388 
1389  // No bounded regions
1390  if (n_regions_bounded_by_facets == 0)
1391  {
1392  if (inner_surface_is_closed)
1393  {
1394  std::ostringstream error_message;
1395  error_message
1396  << "Something fishy going on! "
1397  << "Internal faceted surface " << i_internal
1398  << " is closed but does not bound any regions!\n"
1399  << "Specify one-based region ID for all facets using\n\n"
1400  << " TetMeshFacet::set_one_based_adjacent_region_id(...)\n\n";
1401  throw OomphLibError(error_message.str(),
1402  OOMPH_CURRENT_FUNCTION,
1403  OOMPH_EXCEPTION_LOCATION);
1404  }
1405  }
1406  else
1407  {
1408  // Do we need to create a hole in the main volume that contains
1409  // the multiple sub-volumes/regions?
1410  unsigned offset_for_extra_hole = 0;
1411  if (n_regions_bounded_by_facets > 1)
1412  {
1413  geo_file << std::endl;
1414  geo_file << "// Define Plane Surfaces bounding compound volume"
1415  << std::endl
1416  << "//-----------------------------------------------"
1417  << std::endl;
1418  geo_file << "// that will have to be treated as hole in main volume"
1419  << std::endl
1420  << "//-----------------------------------------------"
1421  << std::endl;
1422  offset_for_extra_hole = 1;
1423  geo_file << "Surface Loop(" << nvolume_offset + 1 << ") = {";
1424  unsigned n = outer_bounding_facet.size();
1425  for (unsigned f = 0; f < n - 1; f++)
1426  {
1427  geo_file << outer_bounding_facet[f] << ",";
1428  }
1429  geo_file << outer_bounding_facet[n - 1] << "};" << std::endl;
1430 
1431  // Bump
1432  number_of_volumes_created_for_this_internal_object++;
1433  }
1434 
1435  // Need to remove the internal volume from the outer volume
1436  volume_id_to_be_subtracted_off.push_back(nvolume_offset + 1);
1437 
1438  // Deal with actual regions that fill the hole
1439  unsigned count = 0;
1440  for (std::map<unsigned, Vector<unsigned>>::iterator it =
1441  region_bounding_facet.begin();
1442  it != region_bounding_facet.end();
1443  it++)
1444  {
1445  geo_file << std::endl;
1446  geo_file << "// Define Plane Surfaces bounding the region volume "
1447  << (*it).first << std::endl;
1448  geo_file << "//----------------------------------------------------"
1449  << std::endl;
1450  geo_file << "Surface Loop("
1451  << nvolume_offset + 1 + offset_for_extra_hole + count
1452  << ") = {";
1453  unsigned n = (*it).second.size();
1454  for (unsigned f = 0; f < n - 1; f++)
1455  {
1456  geo_file << ((*it).second)[f] << ",";
1457  }
1458  geo_file << ((*it).second)[n - 1] << "};" << std::endl;
1459 
1460  geo_file << std::endl;
1461  geo_file << "// Define volume "
1462  << nvolume_offset + 1 + offset_for_extra_hole + count
1463  << " as the volume bounded by Surface Loop "
1464  << nvolume_offset + 1 + offset_for_extra_hole + count
1465  << std::endl;
1466  geo_file
1467  << "//--------------------------------------------------------"
1468  << std::endl;
1469  geo_file << "Volume("
1470  << nvolume_offset + 1 + offset_for_extra_hole + count
1471  << ")={"
1472  << nvolume_offset + 1 + offset_for_extra_hole + count
1473  << "};" << std::endl;
1474  geo_file << std::endl;
1475 
1476  // Bump
1477  number_of_volumes_created_for_this_internal_object++;
1478 
1479 
1480  // Add as physical (to be meshed) volume if it's not a volume
1481  bool mesh_the_volume = true;
1482  if (closed_srf_pt != 0)
1483  {
1484  if (closed_srf_pt->faceted_volume_represents_hole_for_gmsh())
1485  {
1486  mesh_the_volume = false;
1487  }
1488  }
1489  if (mesh_the_volume)
1490  {
1491  geo_file << "// Define one-based region IDs" << std::endl;
1492  geo_file << "//----------------------------" << std::endl;
1493  geo_file << "Physical Volume(" << (*it).first << ")={"
1494  << nvolume_offset + 1 + offset_for_extra_hole + count
1495  << "};" << std::endl;
1496 
1497  unsigned ns_embedded =
1498  surfaces_to_be_embedded_in_specified_one_based_region[(*it)
1499  .first]
1500  .size();
1501  if (ns_embedded > 0)
1502  {
1503  geo_file << "// This region has " << ns_embedded
1504  << " embedded surfaces\n";
1505  geo_file << "Surface{";
1506  for (unsigned i = 0; i < ns_embedded - 1; i++)
1507  {
1508  geo_file
1509  << surfaces_to_be_embedded_in_specified_one_based_region
1510  [(*it).first][i]
1511  << ",";
1512  }
1513  geo_file
1514  << surfaces_to_be_embedded_in_specified_one_based_region
1515  [(*it).first][ns_embedded - 1]
1516  << "}In Volume {"
1517  << nvolume_offset + 1 + offset_for_extra_hole + count << "};"
1518  << std::endl;
1519  }
1520  }
1521 
1522  count++;
1523  }
1524  }
1525 
1526  geo_file << std::endl;
1527  geo_file << "// Define one-based boundary IDs" << std::endl;
1528  geo_file << "//------------------------------" << std::endl;
1529  for (unsigned f = 0; f < nfacet; f++)
1530  {
1531  unsigned one_based_boundary_id =
1532  internal_surface_pt[i_internal]->one_based_facet_boundary_id(f);
1533  geo_file << "Physical Surface(" << one_based_boundary_id << ") = {"
1534  << nfacet_offset + f + 1 << "};" << std::endl;
1535  }
1536 
1537  // Bump
1538  nvertex_offset += internal_surface_pt[i_internal]->nvertex();
1539  nfacet_offset += internal_surface_pt[i_internal]->nfacet();
1540  nvolume_offset += number_of_volumes_created_for_this_internal_object;
1541  nline_offset += tet_edge.size();
1542  }
1543 
1544 
1545  // Done with the internal regions; write the actual volume
1546  // and specify embedded surfaces
1547  geo_file << std::endl;
1548  geo_file << "// Define volume 1 as the volume bounded by Surface Loop 1"
1549  << std::endl;
1550  geo_file << "//--------------------------------------------------------"
1551  << std::endl;
1552  unsigned n = volume_id_to_be_subtracted_off.size();
1553  if (n > 0)
1554  {
1555  geo_file << "// with volume[s]: ";
1556  for (unsigned i = 0; i < n; i++)
1557  {
1558  geo_file << volume_id_to_be_subtracted_off[i] << " ";
1559  }
1560  geo_file << "removed." << std::endl;
1561  geo_file << "//--------------------------------------------------------"
1562  << std::endl;
1563  }
1564 
1565  // Add initial volume
1566  geo_file << "Volume(1)={1";
1567 
1568  // Remove volumes that has separate IDs
1569  for (unsigned i = 0; i < n; i++)
1570  {
1571  geo_file << "," << volume_id_to_be_subtracted_off[i];
1572  }
1573  geo_file << "};" << std::endl;
1574  geo_file << std::endl;
1575 
1576 
1577  geo_file << "// Define one-based region IDs" << std::endl;
1578  geo_file << "//----------------------------" << std::endl;
1579  geo_file << "Physical Volume(1"
1580  << ")={1};" << std::endl;
1581 
1582 
1583  // Now embed any surfaces that don't bound volumes
1584  unsigned ns = surfaces_to_be_embedded_in_main_volume.size();
1585  if (ns > 0)
1586  {
1587  geo_file << std::endl;
1588  geo_file << "// Embed Plane Surfaces in main volume (volume 1)"
1589  << std::endl;
1590  geo_file << "//-----------------------------------------------"
1591  << std::endl;
1592  geo_file << "Surface{";
1593  for (unsigned s = 0; s < ns - 1; s++)
1594  {
1595  geo_file << surfaces_to_be_embedded_in_main_volume[s] << ",";
1596  }
1597  geo_file << surfaces_to_be_embedded_in_main_volume[ns - 1]
1598  << "} In Volume{1};" << std::endl;
1599  geo_file << std::endl;
1600  }
1601 
1602 
1603  // Mesh grading
1604  {
1605  if (use_mesh_grading_from_file)
1606  {
1607 #ifdef PARANOID
1608 
1609  // Open wide...
1610  std::ifstream file(target_size_file_name.c_str(), std::ios_base::in);
1611 
1612  // Check that the file actually opened correctly
1613  if (!file.is_open())
1614  {
1615  std::string error_msg("Failed to open target volume file: ");
1616  error_msg += "\"" + target_size_file_name;
1617  throw OomphLibError(
1618  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1619  }
1620 #endif
1621 
1622  geo_file << "Field[1]=Structured;" << std::endl;
1623  geo_file << "Field[1].FileName=\"" << target_size_file_name << "\";"
1624  << std::endl;
1625  geo_file << "Field[1].TextFormat=1;" << std::endl;
1626  geo_file << "Background Field = 1;" << std::endl;
1627  }
1628  }
1629 
1630  // Mesh the bloody thing
1631  geo_file << std::endl;
1632  geo_file << "Mesh 3;" << std::endl;
1633 
1634 
1635  /* // hierher Christmas attempt at optimisation */
1636  /* geo_file << std::endl; */
1637  /* geo_file << "// Attempt at optimisation:
1638  * http://onelab.info/pipermail/gmsh/2015/010126.html" << std::endl; */
1639  /* geo_file << "//----------------------------" << std::endl; */
1640  /* geo_file << "Mesh.Optimize=1;" << std::endl; */
1641  /* geo_file << "Mesh.OptimizeNetgen=1;" << std::endl; */
1642 
1643 
1644  geo_file.close();
1645  }
1646 
1647  /// Parameters
1649  };
1650 
1651 
1652  /// ///////////////////////////////////////////////////////////////////////
1653  /// ///////////////////////////////////////////////////////////////////////
1654  /// ///////////////////////////////////////////////////////////////////////
1655 
1656 
1657  //=========================================================================
1658  // Gmsh-based Tet Mesh
1659  //=========================================================================
1660  template<class ELEMENT>
1661  class GmshTetMesh : public virtual TetMeshBase, public virtual Mesh
1662  {
1663  public:
1664  /// Constructor
1665  GmshTetMesh(GmshParameters* gmsh_parameters_pt,
1666  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1667  : Gmsh_parameters_pt(gmsh_parameters_pt)
1668  {
1669  bool use_mesh_grading_from_file = false;
1670  build_it(time_stepper_pt, use_mesh_grading_from_file);
1671  }
1672 
1673  /// Constructor. If boolean is set
1674  /// to true, the target element sizes are read from file (used during
1675  /// adaptation; otherwise uniform target size is used).
1676  GmshTetMesh(GmshParameters* gmsh_parameters_pt,
1677  const bool& use_mesh_grading_from_file,
1678  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1679  : Gmsh_parameters_pt(gmsh_parameters_pt)
1680  {
1681  build_it(time_stepper_pt, use_mesh_grading_from_file);
1682  }
1683 
1684 
1685  private:
1686  // Build function
1687  void build_it(TimeStepper* time_stepper_pt,
1688  const bool& use_mesh_grading_from_file)
1689  {
1690  // Mesh can only be built with 3D Telements.
1691  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
1692 
1693  // Transfer data from parameters
1694  TetMeshFacetedClosedSurface* outer_boundary_pt =
1696 
1697  Vector<TetMeshFacetedSurface*> internal_surface_pt =
1699 
1700  // Remember timestepper
1701  Time_stepper_pt = time_stepper_pt;
1702 
1703  // Store the boundary
1704  Outer_boundary_pt = outer_boundary_pt;
1705 
1706  // Setup reverse lookup scheme
1707  {
1708  unsigned n_facet = Outer_boundary_pt->nfacet();
1709  for (unsigned f = 0; f < n_facet; f++)
1710  {
1711  unsigned b = Outer_boundary_pt->one_based_facet_boundary_id(f);
1712  if (b != 0)
1713  {
1714  Tet_mesh_faceted_surface_pt[b - 1] = Outer_boundary_pt;
1715  Tet_mesh_facet_pt[b - 1] = Outer_boundary_pt->facet_pt(f);
1716  }
1717  else
1718  {
1719  std::ostringstream error_message;
1720  error_message << "Boundary IDs have to be one-based. Yours is " << b
1721  << "\n";
1722  throw OomphLibError(error_message.str(),
1723  OOMPH_CURRENT_FUNCTION,
1724  OOMPH_EXCEPTION_LOCATION);
1725  }
1726  }
1727  }
1728 
1729 
1730  // Store the internal boundary
1731  Internal_surface_pt = internal_surface_pt;
1732 
1733  // Setup reverse lookup scheme
1734  {
1735  unsigned n = Internal_surface_pt.size();
1736  for (unsigned i = 0; i < n; i++)
1737  {
1738  unsigned n_facet = Internal_surface_pt[i]->nfacet();
1739  for (unsigned f = 0; f < n_facet; f++)
1740  {
1741  unsigned b = Internal_surface_pt[i]->one_based_facet_boundary_id(f);
1742  if (b != 0)
1743  {
1744  Tet_mesh_faceted_surface_pt[b - 1] = Internal_surface_pt[i];
1745  Tet_mesh_facet_pt[b - 1] = Internal_surface_pt[i]->facet_pt(f);
1746  }
1747  else
1748  {
1749  std::ostringstream error_message;
1750  error_message << "Boundary IDs have to be one-based. Yours is "
1751  << b << "\n";
1752  throw OomphLibError(error_message.str(),
1753  OOMPH_CURRENT_FUNCTION,
1754  OOMPH_EXCEPTION_LOCATION);
1755  }
1756  }
1757  }
1758  }
1759 
1760  // Build scaffold
1761  GmshTetScaffoldMesh* tmp_scaffold_mesh_pt =
1762  new GmshTetScaffoldMesh(Gmsh_parameters_pt, use_mesh_grading_from_file);
1763 
1764  // Convert mesh from scaffold to actual mesh
1765  build_from_scaffold(tmp_scaffold_mesh_pt, time_stepper_pt);
1766 
1767  // Kill the scaffold
1768  delete tmp_scaffold_mesh_pt;
1769  tmp_scaffold_mesh_pt = 0;
1770 
1771  // Setup boundary coordinates
1772  unsigned nb = nboundary();
1773  for (unsigned b = 0; b < nb; b++)
1774  {
1775  bool switch_normal = false;
1776  setup_boundary_coordinates<ELEMENT>(b, switch_normal);
1777  }
1778 
1779  // Now snap onto geometric objects associated with triangular facets
1780  // (if any!)
1781  snap_nodes_onto_geometric_objects();
1782  }
1783 
1784 
1785  protected:
1786  /// Parameters
1788 
1789 
1790  private:
1791  /// Build unstructured tet gmesh mesh based on output from scaffold
1792  void build_from_scaffold(GmshTetScaffoldMesh* tmp_scaffold_mesh_pt,
1793  TimeStepper* time_stepper_pt)
1794  {
1795  // Mesh can only be built with 3D Telements.
1796  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
1797 
1798  // Create space for elements
1799  unsigned nelem = tmp_scaffold_mesh_pt->nelement();
1800  Element_pt.resize(nelem);
1801 
1802  // Relation between elements pointers and numbers in old mesh
1803  std::map<FiniteElement*, unsigned> scaffold_mesh_element_number;
1804 
1805  // Create space for nodes
1806  unsigned nnode_scaffold = tmp_scaffold_mesh_pt->nnode();
1807  Node_pt.resize(nnode_scaffold);
1808 
1809  // Set number of boundaries
1810  unsigned nbound = tmp_scaffold_mesh_pt->nboundary();
1811  set_nboundary(nbound);
1812 
1813  // Resize boundary info stuff
1814  Boundary_element_pt.resize(nbound); // hierher shouldn't this be a map?
1815  Face_index_at_boundary.resize(nbound); // hierher shouldn't this be a map?
1816  Boundary_region_element_pt.resize(nbound);
1817  Face_index_region_at_boundary.resize(nbound);
1818 
1819  // Build elements
1820  for (unsigned e = 0; e < nelem; e++)
1821  {
1822  Element_pt[e] = new ELEMENT;
1823  }
1824 
1825  // Number of nodes per element
1826  unsigned nnod_el = tmp_scaffold_mesh_pt->finite_element_pt(0)->nnode();
1827 
1828  // Setup map to check the (pseudo-)global node number
1829  // Nodes whose number is zero haven't been copied across
1830  // into the mesh yet.
1831  std::map<Node*, unsigned> global_number;
1832  unsigned global_count = 0;
1833 
1834  // Loop over elements in scaffold mesh, visit their nodes
1835  for (unsigned e = 0; e < nelem; e++)
1836  {
1837  // Setup reverse lookup scheme to decipher lookup schemes from
1838  // scaffold mesh
1839  scaffold_mesh_element_number[tmp_scaffold_mesh_pt->finite_element_pt(
1840  e)] = e;
1841 
1842  // Loop over all nodes in element
1843  for (unsigned j = 0; j < nnod_el; j++)
1844  {
1845  // Pointer to node in the scaffold mesh
1846  Node* scaffold_node_pt =
1847  tmp_scaffold_mesh_pt->finite_element_pt(e)->node_pt(j);
1848 
1849  // Get the (pseudo-)global node number in scaffold mesh
1850  // (It's zero [=default] if not visited this one yet)
1851  unsigned j_global = global_number[scaffold_node_pt];
1852 
1853  // Haven't done this one yet
1854  if (j_global == 0)
1855  {
1856  // Get pointer to set of mesh boundaries that this
1857  // scaffold node occupies; NULL if the node is not on any boundary
1858  std::set<unsigned>* boundaries_pt;
1859  scaffold_node_pt->get_boundaries_pt(boundaries_pt);
1860 
1861  // Is it on boundaries?
1862  if (boundaries_pt != 0)
1863  {
1864  // Create new boundary node
1865  Node* new_node_pt = finite_element_pt(e)->construct_boundary_node(
1866  j, time_stepper_pt);
1867 
1868  // Give it a number (not necessarily the global node
1869  // number in the scaffold mesh -- we just need something
1870  // to keep track...)
1871  global_count++;
1872  global_number[scaffold_node_pt] = global_count;
1873 
1874  // Add to boundaries
1875  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
1876  it != boundaries_pt->end();
1877  ++it)
1878  {
1879  add_boundary_node(*it, new_node_pt);
1880  }
1881  }
1882  // Build normal node
1883  else
1884  {
1885  // Create new normal node
1886  finite_element_pt(e)->construct_node(j, time_stepper_pt);
1887 
1888  // Give it a number (not necessarily the global node
1889  // number in the scaffold mesh -- we just need something
1890  // to keep track...)
1891  global_count++;
1892  global_number[scaffold_node_pt] = global_count;
1893  }
1894 
1895  // Copy new node, created using the NEW element's construct_node
1896  // function into global storage, using the same global
1897  // node number that we've just associated with the
1898  // corresponding node in the scaffold mesh
1899  Node_pt[global_count - 1] = finite_element_pt(e)->node_pt(j);
1900 
1901  // Assign coordinates
1902  Node_pt[global_count - 1]->x(0) = scaffold_node_pt->x(0);
1903  Node_pt[global_count - 1]->x(1) = scaffold_node_pt->x(1);
1904  Node_pt[global_count - 1]->x(2) = scaffold_node_pt->x(2);
1905  }
1906  // This one has already been done: Copy across
1907  else
1908  {
1909  finite_element_pt(e)->node_pt(j) = Node_pt[j_global - 1];
1910  }
1911  }
1912 
1913 
1914  // Now figure out which boundaries the faces are on
1915  FiniteElement* fe_pt = finite_element_pt(e);
1916  for (unsigned f = 0; f < 4; f++)
1917  {
1918  Node* face_node0_pt = 0;
1919  Node* face_node1_pt = 0;
1920  Node* face_node2_pt = 0;
1921 
1922  switch (f)
1923  {
1924  case 0:
1925  face_node0_pt = fe_pt->node_pt(1);
1926  face_node1_pt = fe_pt->node_pt(2);
1927  face_node2_pt = fe_pt->node_pt(3);
1928  break;
1929 
1930  case 1:
1931  face_node0_pt = fe_pt->node_pt(0);
1932  face_node1_pt = fe_pt->node_pt(2);
1933  face_node2_pt = fe_pt->node_pt(3);
1934  break;
1935 
1936  case 2:
1937  face_node0_pt = fe_pt->node_pt(0);
1938  face_node1_pt = fe_pt->node_pt(1);
1939  face_node2_pt = fe_pt->node_pt(3);
1940  break;
1941 
1942  case 3:
1943  face_node0_pt = fe_pt->node_pt(0);
1944  face_node1_pt = fe_pt->node_pt(1);
1945  face_node2_pt = fe_pt->node_pt(2);
1946  break;
1947 
1948  default:
1949 
1950  std::ostringstream error_message;
1951  error_message << "Wrong face number " << f << std::endl;
1952  throw OomphLibError(error_message.str(),
1953  OOMPH_CURRENT_FUNCTION,
1954  OOMPH_EXCEPTION_LOCATION);
1955  }
1956 
1957  // If any of the boundary sets are empty we don't have
1958  // three nodes on the boundary...
1959  std::set<unsigned>* bc0_pt;
1960  face_node0_pt->get_boundaries_pt(bc0_pt);
1961  if (bc0_pt != 0)
1962  {
1963  std::set<unsigned>* bc1_pt;
1964  face_node1_pt->get_boundaries_pt(bc1_pt);
1965  if (bc1_pt != 0)
1966  {
1967  std::set<unsigned>* bc2_pt;
1968  face_node2_pt->get_boundaries_pt(bc2_pt);
1969  if (bc2_pt != 0)
1970  {
1971  std::set<unsigned> common_bound_0_and_1;
1972  std::set_intersection(
1973  bc0_pt->begin(),
1974  bc0_pt->end(),
1975  bc1_pt->begin(),
1976  bc1_pt->end(),
1977  std::inserter(common_bound_0_and_1,
1978  common_bound_0_and_1.begin()));
1979  std::set<unsigned> common_bound_0_and_1_and_2;
1980  std::set_intersection(
1981  common_bound_0_and_1.begin(),
1982  common_bound_0_and_1.end(),
1983  bc2_pt->begin(),
1984  bc2_pt->end(),
1985  std::inserter(common_bound_0_and_1_and_2,
1986  common_bound_0_and_1_and_2.begin()));
1987  for (std::set<unsigned>::iterator it =
1988  common_bound_0_and_1_and_2.begin();
1989  it != common_bound_0_and_1_and_2.end();
1990  it++)
1991  {
1992  Boundary_element_pt[(*it)].push_back(fe_pt);
1993  Face_index_at_boundary[(*it)].push_back(f);
1994  }
1995  }
1996  }
1997  }
1998  }
1999  }
2000 
2001  // Copy across region information (scaffold mesh is a friend)
2002  unsigned nr = tmp_scaffold_mesh_pt->Region_element_pt.size();
2003  Region_attribute.resize(nr);
2004  Region_element_pt.resize(nr);
2005  for (unsigned i = 0; i < nr; i++)
2006  { //--
2007  Region_attribute[i] = tmp_scaffold_mesh_pt->Region_attribute[i];
2008  unsigned nel = tmp_scaffold_mesh_pt->Region_element_pt[i].size();
2009  Region_element_pt[i].resize(nel);
2010  for (unsigned e = 0; e < nel; e++)
2011  {
2012  FiniteElement* scaff_el_pt =
2013  tmp_scaffold_mesh_pt->Region_element_pt[i][e];
2014  unsigned scaff_el_number = scaffold_mesh_element_number[scaff_el_pt];
2015  Region_element_pt[i][e] =
2016  dynamic_cast<FiniteElement*>(Element_pt[scaff_el_number]);
2017 
2018 
2019  // Now figure out which boundaries the faces are on (again!)
2020  FiniteElement* fe_pt = Region_element_pt[i][e];
2021  for (unsigned f = 0; f < 4; f++)
2022  {
2023  Node* face_node0_pt = 0;
2024  Node* face_node1_pt = 0;
2025  Node* face_node2_pt = 0;
2026 
2027  switch (f)
2028  {
2029  case 0:
2030  face_node0_pt = fe_pt->node_pt(1);
2031  face_node1_pt = fe_pt->node_pt(2);
2032  face_node2_pt = fe_pt->node_pt(3);
2033  break;
2034 
2035  case 1:
2036  face_node0_pt = fe_pt->node_pt(0);
2037  face_node1_pt = fe_pt->node_pt(2);
2038  face_node2_pt = fe_pt->node_pt(3);
2039  break;
2040 
2041  case 2:
2042  face_node0_pt = fe_pt->node_pt(0);
2043  face_node1_pt = fe_pt->node_pt(1);
2044  face_node2_pt = fe_pt->node_pt(3);
2045  break;
2046 
2047  case 3:
2048  face_node0_pt = fe_pt->node_pt(0);
2049  face_node1_pt = fe_pt->node_pt(1);
2050  face_node2_pt = fe_pt->node_pt(2);
2051  break;
2052 
2053  default:
2054  std::ostringstream error_message;
2055  error_message << "Wrong face number " << f << std::endl;
2056  throw OomphLibError(error_message.str(),
2057  OOMPH_CURRENT_FUNCTION,
2058  OOMPH_EXCEPTION_LOCATION);
2059  }
2060 
2061  // If any of the boundary sets are empty we don't have
2062  // three nodes on the boundary...
2063  std::set<unsigned>* bc0_pt;
2064  face_node0_pt->get_boundaries_pt(bc0_pt);
2065  if (bc0_pt != 0)
2066  {
2067  std::set<unsigned>* bc1_pt;
2068  face_node1_pt->get_boundaries_pt(bc1_pt);
2069  if (bc1_pt != 0)
2070  {
2071  std::set<unsigned>* bc2_pt;
2072  face_node2_pt->get_boundaries_pt(bc2_pt);
2073  if (bc2_pt != 0)
2074  {
2075  std::set<unsigned> common_bound_0_and_1;
2076  std::set_intersection(
2077  bc0_pt->begin(),
2078  bc0_pt->end(),
2079  bc1_pt->begin(),
2080  bc1_pt->end(),
2081  std::inserter(common_bound_0_and_1,
2082  common_bound_0_and_1.begin()));
2083  std::set<unsigned> common_bound_0_and_1_and_2;
2084  std::set_intersection(
2085  common_bound_0_and_1.begin(),
2086  common_bound_0_and_1.end(),
2087  bc2_pt->begin(),
2088  bc2_pt->end(),
2089  std::inserter(common_bound_0_and_1_and_2,
2090  common_bound_0_and_1_and_2.begin()));
2091  for (std::set<unsigned>::iterator it =
2092  common_bound_0_and_1_and_2.begin();
2093  it != common_bound_0_and_1_and_2.end();
2094  it++)
2095  {
2096  Boundary_region_element_pt[(*it)][Region_attribute[i]]
2097  .push_back(fe_pt);
2098  Face_index_region_at_boundary[(*it)][Region_attribute[i]]
2099  .push_back(f);
2100  }
2101  }
2102  }
2103  }
2104  }
2105  }
2106  }
2107 
2108  // Lookup scheme has now been setup
2109  Lookup_for_elements_next_boundary_is_setup = true;
2110 
2111 
2112  // At this point we've created all the elements and
2113  // created their vertex nodes. Now we need to create
2114  // the additional (midside and internal) nodes!
2115 
2116  // Get number of nodes along element edge
2117  unsigned n_node_1d = finite_element_pt(0)->nnode_1d();
2118 
2119  // At the moment we're only able to deal with nnode_1d=2 or 3.
2120  if ((n_node_1d != 2) && (n_node_1d != 3))
2121  {
2122  std::ostringstream error_message;
2123  error_message << "Mesh generation from gmsh currently only works\n";
2124  error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
2125  error_message << "for nnode_1d=" << n_node_1d << std::endl;
2126  throw OomphLibError(error_message.str(),
2127  OOMPH_CURRENT_FUNCTION,
2128  OOMPH_EXCEPTION_LOCATION);
2129  }
2130 
2131  // Storage for the local coordinate of the new node
2132  Vector<double> s(3);
2133 
2134  // Map associating edge with midside node
2135  std::map<std::pair<Node*, Node*>, Node*> midside_node_pt;
2136 
2137  // Loop over all elements
2138  for (unsigned e = 0; e < nelem; e++)
2139  {
2140  // Cache pointers to the elements
2141  FiniteElement* const el_pt = this->finite_element_pt(e);
2142  FiniteElement* const simplex_el_pt =
2143  tmp_scaffold_mesh_pt->finite_element_pt(e);
2144 
2145  // Loop over the edges
2146  for (unsigned j = 0; j < 6; ++j)
2147  {
2148  Node* nod0_pt = 0;
2149  Node* nod1_pt = 0;
2150  unsigned new_node_number = 0;
2151  std::pair<Node*, Node*> edge;
2152  switch (j)
2153  {
2154  case 0:
2155  nod0_pt = el_pt->node_pt(0);
2156  nod1_pt = el_pt->node_pt(1);
2157  new_node_number = 4;
2158  break;
2159 
2160  case 1:
2161  nod0_pt = el_pt->node_pt(0);
2162  nod1_pt = el_pt->node_pt(2);
2163  new_node_number = 5;
2164  break;
2165 
2166  case 2:
2167  nod0_pt = el_pt->node_pt(0);
2168  nod1_pt = el_pt->node_pt(3);
2169  new_node_number = 6;
2170  break;
2171 
2172  case 3:
2173  nod0_pt = el_pt->node_pt(1);
2174  nod1_pt = el_pt->node_pt(2);
2175  new_node_number = 7;
2176  break;
2177 
2178  case 4:
2179  nod0_pt = el_pt->node_pt(2);
2180  nod1_pt = el_pt->node_pt(3);
2181  new_node_number = 8;
2182  break;
2183 
2184  case 5:
2185  nod0_pt = el_pt->node_pt(1);
2186  nod1_pt = el_pt->node_pt(3);
2187  new_node_number = 9;
2188  break;
2189 
2190  default:
2191  std::ostringstream error_message;
2192  error_message << "Wrong edge number " << j << std::endl;
2193  throw OomphLibError(error_message.str(),
2194  OOMPH_CURRENT_FUNCTION,
2195  OOMPH_EXCEPTION_LOCATION);
2196  }
2197 
2198  // Identify existence of node via edge
2199  edge = std::make_pair(std::min(nod0_pt, nod1_pt),
2200  std::max(nod0_pt, nod1_pt));
2201  Node* existing_node_pt = midside_node_pt[edge];
2202  if (existing_node_pt == 0)
2203  {
2204  // If any of the boundary sets are empty we don't have
2205  // two edge nodes on the boundary...
2206  std::set<unsigned> common_bound_0_and_1;
2207  std::set<unsigned>* bc0_pt;
2208  nod0_pt->get_boundaries_pt(bc0_pt);
2209  if (bc0_pt != 0)
2210  {
2211  std::set<unsigned>* bc1_pt;
2212  nod1_pt->get_boundaries_pt(bc1_pt);
2213  if (bc1_pt != 0)
2214  {
2215  std::set_intersection(
2216  bc0_pt->begin(),
2217  bc0_pt->end(),
2218  bc1_pt->begin(),
2219  bc1_pt->end(),
2220  std::inserter(common_bound_0_and_1,
2221  common_bound_0_and_1.begin()));
2222  }
2223  }
2224  // Storage for the new node
2225  Node* new_node_pt = 0;
2226 
2227  // New non-boundary node:
2228  if (common_bound_0_and_1.size() == 0)
2229  {
2230  new_node_pt =
2231  el_pt->construct_node(new_node_number, time_stepper_pt);
2232  }
2233  // New boundary node
2234  else
2235  {
2236  new_node_pt = el_pt->construct_boundary_node(new_node_number,
2237  time_stepper_pt);
2238  for (std::set<unsigned>::iterator it =
2239  common_bound_0_and_1.begin();
2240  it != common_bound_0_and_1.end();
2241  it++)
2242  {
2243  this->add_boundary_node((*it), new_node_pt);
2244  }
2245  }
2246 
2247  // Find the local coordinates of the node
2248  el_pt->local_coordinate_of_node(new_node_number, s);
2249 
2250  // Find the coordinates of the new node from the existing
2251  // and fully-functional element in the scaffold mesh
2252  for (unsigned i = 0; i < 3; i++)
2253  {
2254  new_node_pt->x(i) = simplex_el_pt->interpolated_x(s, i);
2255  }
2256 
2257  // Associate node with edge
2258  midside_node_pt[edge] = new_node_pt;
2259 
2260  // Add node to mesh
2261  Node_pt.push_back(new_node_pt);
2262  }
2263  // Node already exists
2264  else
2265  {
2266  el_pt->node_pt(new_node_number) = existing_node_pt;
2267  }
2268  }
2269  }
2270  }
2271  };
2272 
2273 
2274  /// ///////////////////////////////////////////////////////////////////////
2275  /// ///////////////////////////////////////////////////////////////////////
2276  /// ///////////////////////////////////////////////////////////////////////
2277 
2278 
2279  //=========================================================================
2280  // Refineable Gmsh-based Tet Mesh
2281  //=========================================================================
2282  template<class ELEMENT>
2283  class RefineableGmshTetMesh : public virtual GmshTetMesh<ELEMENT>,
2284  public virtual RefineableTetMeshBase
2285  {
2286  public:
2287  /// Constructor. If boolean is set
2288  /// to true, the target element sizes are read from file (used during
2289  /// adaptation; otherwise uniform target size is used).
2291  GmshParameters* gmsh_parameters_pt,
2292  const bool& use_mesh_grading_from_file,
2293  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2294  : GmshTetMesh<ELEMENT>(
2295  gmsh_parameters_pt, use_mesh_grading_from_file, time_stepper_pt)
2296  {
2298  }
2299 
2300  /// Constructor
2302  GmshParameters* gmsh_parameters_pt,
2303  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2304  : GmshTetMesh<ELEMENT>(gmsh_parameters_pt,
2305  false, // Don't read mesh size data from file
2306  // unless explicitly requested.
2307  time_stepper_pt)
2308  {
2310  }
2311 
2312 
2313  /// Adapt mesh, based on elemental error provided
2314  void adapt(const Vector<double>& elem_error);
2315 
2316  /// Refine uniformly
2318  {
2319  // hierher do it!
2320  std::string error_msg("Not written yet...");
2321  throw OomphLibError(
2322  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2323  }
2324 
2325  /// Unrefine uniformly
2326  void refine_uniformly(DocInfo& doc_info)
2327  {
2328  // hierher do it!
2329  std::string error_msg("Not written yet...");
2330  throw OomphLibError(
2331  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2332  }
2333 
2334 
2335  protected:
2336  /// Helper function to initialise data associated with adaptation
2338  {
2339  // Set max and min targets for adaptation
2340  this->Max_element_size = this->Gmsh_parameters_pt->max_element_size();
2341  this->Min_element_size = this->Gmsh_parameters_pt->min_element_size();
2342  this->Max_permitted_edge_ratio =
2344  }
2345  };
2346 
2347 
2348  /// //////////////////////////////////////////////////////////////////////////
2349  /// //////////////////////////////////////////////////////////////////////////
2350  /// //////////////////////////////////////////////////////////////////////////
2351 
2352 
2353 } // namespace oomph
2354 
2355 #endif
Class to collate parameters for Gmsh mesh generation.
double Dx_for_target_size_transfer
(Isotropic) grid spacing for target size transfer
Vector< TetMeshFacetedSurface * > & internal_surface_pt()
Internal boundaries.
TetMeshFacetedClosedSurface * Outer_boundary_pt
Pointer to outer boundary.
std::string Geo_and_msh_file_stem
Stem for geo and msh files (input/output to command-line gmsh invocation)
double & min_element_size()
Min. element size during refinement.
double & element_volume()
Uniform target element volume.
std::string Gmsh_command_line_invocation
Gmsh command line invocation string.
std::string & stem_for_filename_gmsh_size_transfer()
Stem for filename used to doc target element sizes on gmsh grid. No doc if stem is equal to empty str...
double Max_element_size
Max. element size during refinement.
GmshParameters(TetMeshFacetedClosedSurface *const &outer_boundary_pt, const std::string &gmsh_command_line_invocation)
Specify outer boundary of domain to be meshed. Other parameters get default values and can be set via...
unsigned Max_sample_points_for_limited_locate_zeta_during_target_size_transfer
Target size is transferred onto regular grid (for gmsh) by locate zeta. We try to find the exact poin...
unsigned & max_sample_points_for_limited_locate_zeta_during_target_size_transfer()
Target size is transferred onto regular grid (for gmsh) by locate zeta. We try to find the exact poin...
std::string & target_size_file_name()
Filename for target volumes (for system-call based transfer to gmsh during mesh adaptation)....
double Min_element_size
Min. element size during refinement.
double Element_volume
Uniform element volume.
Vector< TetMeshFacetedSurface * > Internal_surface_pt
Internal boundaries.
std::string & gmsh_command_line_invocation()
String to be issued via system command to activate gmsh.
unsigned Gmsh_onscreen_output_counter
Counter for marker that indicates where we are in gmsh on-screen output.
std::string & geo_and_msh_file_stem()
Stem for geo and msh files (input/output to command-line gmsh invocation)
TetMeshFacetedClosedSurface *& outer_boundary_pt()
Outer boundary.
double Max_permitted_edge_ratio
Max edge ratio before remesh gets triggered.
std::string & gmsh_onscreen_output_file_name()
Output filename for gmsh on-screen output.
bool projection_is_disabled()
Is projection of old solution onto new mesh disabled?
int & counter_for_filename_gmsh_size_transfer()
Counter for filename used to doc target element sizes on gmsh grid. No doc if stem is equal to empty ...
std::string Stem_for_filename_gmsh_size_transfer
Stem for filename used to doc target element sizes on gmsh grid. No doc if stem is equal to empty str...
void disable_projection()
Disable projection of old solution onto new mesh.
std::string Target_size_file_name
Filename for target volume (for system-call based transfer to gmsh during mesh adaptation)
double & dx_for_target_size_transfer()
(Isotropic) grid spacing for target size transfer
bool Projection_is_disabled
Is projection of old solution onto new mesh disabled?
std::string Gmsh_onscreen_output_file_name
Output filename for gmsh on-screen output.
int Counter_for_filename_gmsh_size_transfer
Counter for filename used to doc target element sizes on gmsh grid. No doc if stem is equal to empty ...
void enable_projection()
Disable projection of old solution onto new mesh.
double & max_element_size()
Max. element size during refinement.
unsigned & gmsh_onscreen_output_counter()
Counter for marker that indicates where we are in gmsh on-screen output.
double & max_permitted_edge_ratio()
Max. permitted edge ratio.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
GmshTetMesh(GmshParameters *gmsh_parameters_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor.
GmshParameters * Gmsh_parameters_pt
Parameters.
void build_from_scaffold(GmshTetScaffoldMesh *tmp_scaffold_mesh_pt, TimeStepper *time_stepper_pt)
Build unstructured tet gmesh mesh based on output from scaffold.
GmshTetMesh(GmshParameters *gmsh_parameters_pt, const bool &use_mesh_grading_from_file, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor. If boolean is set to true, the target element sizes are read from file (used during adap...
void build_it(TimeStepper *time_stepper_pt, const bool &use_mesh_grading_from_file)
void create_mesh_from_msh_file()
Create mesh from msh file (created internally via disk-based operations)
GmshTetScaffoldMesh(GmshParameters *gmsh_parameters_pt, const bool &use_mesh_grading_from_file)
Build mesh, based on specified parameters. If boolean is set to true, the target element sizes are re...
void write_geo_file(const bool &use_mesh_grading_from_file)
Write geo file for gmsh.
GmshParameters * Gmsh_parameters_pt
Parameters.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
void refine_uniformly(DocInfo &doc_info)
Unrefine uniformly.
unsigned unrefine_uniformly()
Refine uniformly.
RefineableGmshTetMesh(GmshParameters *gmsh_parameters_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor.
void adapt(const Vector< double > &elem_error)
Adapt mesh, based on elemental error provided.
RefineableGmshTetMesh(GmshParameters *gmsh_parameters_pt, const bool &use_mesh_grading_from_file, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor. If boolean is set to true, the target element sizes are read from file (used during adap...
void initialise_adaptation_data()
Helper function to initialise data associated with adaptation.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
bool operator==(const TetEdge &tet_edge) const
Comparison operator: Edges are identical if their sorted (and therefore possibly reversed) vertex ids...
std::pair< unsigned, unsigned > Vertex_pair
The vertices (sorted by vertex ids)
unsigned second_vertex_id() const
Second vertex id.
bool is_reversed() const
Edge is reversed in the sense that vertex1 actually has a higher id than vertex2 (when specified in t...
TetEdge(const unsigned &vertex1, const unsigned &vertex2)
Constructor: Pass two vertices, identified by their indices Edge "direction" is from lower vertex to ...
bool operator<(const TetEdge &tet_edge) const
Comparison operator. Lexicographic comparison based on vertex ids.
unsigned first_vertex_id() const
First vertex id.
bool Reversed
Is it reversed? I.e. is the first input vertex stored after the second one?
////////////////////////////////////////////////////////////////////// //////////////////////////////...