Public Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
oomph::GmshTetScaffoldMesh Class Reference

#include <gmsh_tet_mesh.template.h>

+ Inheritance diagram for oomph::GmshTetScaffoldMesh:

Public Member Functions

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

Private Member Functions

void create_mesh_from_msh_file ()
 Create mesh from msh file (created internally via disk-based operations) More...
 
void write_geo_file (const bool &use_mesh_grading_from_file)
 Write geo file for gmsh. More...
 

Private Attributes

GmshParameters * Gmsh_parameters_pt
 Parameters. More...
 

Friends

template<class ELEMENT >
class GmshTetMesh
 We're friends with the actual mesh. More...
 

Additional Inherited Members

- Public Types inherited from oomph::Mesh
typedef void(FiniteElement::* SteadyExactSolutionFctPt) (const Vector< double > &x, Vector< double > &soln)
 
typedef void(FiniteElement::* UnsteadyExactSolutionFctPt) (const double &time, const Vector< double > &x, Vector< double > &soln)
 
- Static Public Attributes inherited from oomph::TetMeshBase
static double Tolerance_for_boundary_finding = 1.0e-5
 
- Static Public Attributes inherited from oomph::Mesh
static Steady< 0 > Default_TimeStepper
 The Steady Timestepper. More...
 
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
 Static boolean flag to control warning about mesh level timesteppers. More...
 
- Protected Member Functions inherited from oomph::Mesh
unsigned long assign_global_eqn_numbers (Vector< double * > &Dof_pt)
 Assign (global) equation numbers to the nodes. More...
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign local equation numbers in all elements. More...
 
void convert_to_boundary_node (Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
 
void convert_to_boundary_node (Node *&node_pt)
 
- Protected Attributes inherited from oomph::TetMeshBase
Vector< Vector< FiniteElement * > > Region_element_pt
 
Vector< double > Region_attribute
 
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
 Storage for elements adjacent to a boundary in a particular region. More...
 
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
 
TetMeshFacetedClosedSurface * Outer_boundary_pt
 Faceted surface that defines outer boundaries. More...
 
Vector< TetMeshFacetedSurface * > Internal_surface_pt
 Vector to faceted surfaces that define internal boundaries. More...
 
std::map< unsigned, TetMeshFacetedSurface * > Tet_mesh_faceted_surface_pt
 
std::map< unsigned, TetMeshFacet * > Tet_mesh_facet_pt
 
std::map< unsigned, Vector< Vector< double > > > Triangular_facet_vertex_boundary_coordinate
 
TimeStepper * Time_stepper_pt
 Timestepper used to build nodes. More...
 
- Protected Attributes inherited from oomph::Mesh
Vector< Vector< Node * > > Boundary_node_pt
 
bool Lookup_for_elements_next_boundary_is_setup
 
Vector< Vector< FiniteElement * > > Boundary_element_pt
 
Vector< Vector< int > > Face_index_at_boundary
 
Vector< Node * > Node_pt
 Vector of pointers to nodes. More...
 
Vector< GeneralisedElement * > Element_pt
 Vector of pointers to generalised elements. More...
 
std::vector< bool > Boundary_coordinate_exists
 

Constructor & Destructor Documentation

◆ GmshTetScaffoldMesh()

oomph::GmshTetScaffoldMesh::GmshTetScaffoldMesh ( GmshParameters *  gmsh_parameters_pt,
const bool &  use_mesh_grading_from_file 
)
inline

Build mesh, based on specified parameters. If boolean is set to true, the target element sizes are read from file (used during adaptation; otherwise uniform target size is used).

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 +=
414  Gmsh_parameters_pt->gmsh_command_line_invocation() + " " +
415  Gmsh_parameters_pt->geo_and_msh_file_stem() + ".geo -3";
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
432  create_mesh_from_msh_file();
433 
434  oomph_info << "Time for creating mesh from msh file: "
435  << TimingHelpers::timer() - t_start << " sec " << std::endl;
436  }
oomph::GmshParameters::gmsh_command_line_invocation
std::string & gmsh_command_line_invocation()
String to be issued via system command to activate gmsh.
Definition: gmsh_tet_mesh.template.h:103
oomph::GmshParameters::geo_and_msh_file_stem
std::string & geo_and_msh_file_stem()
Definition: gmsh_tet_mesh.template.h:110
oomph::GmshTetScaffoldMesh::create_mesh_from_msh_file
void create_mesh_from_msh_file()
Create mesh from msh file (created internally via disk-based operations)
Definition: gmsh_tet_mesh.template.h:440
oomph::GmshTetScaffoldMesh::write_geo_file
void write_geo_file(const bool &use_mesh_grading_from_file)
Write geo file for gmsh.
Definition: gmsh_tet_mesh.template.h:1020
oomph::GmshTetScaffoldMesh::Gmsh_parameters_pt
GmshParameters * Gmsh_parameters_pt
Parameters.
Definition: gmsh_tet_mesh.template.h:1648
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
oomph::TimingHelpers::timer
double timer()
returns the time in seconds after some point in past
Definition: oomph_utilities.cc:1295
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References create_mesh_from_msh_file(), oomph::GmshParameters::geo_and_msh_file_stem(), oomph::GmshParameters::gmsh_command_line_invocation(), oomph::GmshParameters::gmsh_onscreen_output_counter(), oomph::GmshParameters::gmsh_onscreen_output_file_name(), Gmsh_parameters_pt, oomph::oomph_info, oomph::Global_string_for_annotation::string(), oomph::TimingHelpers::timer(), and write_geo_file().

Member Function Documentation

◆ create_mesh_from_msh_file()

void oomph::GmshTetScaffoldMesh::create_mesh_from_msh_file ( )
inlineprivate

Create mesh from msh file (created internally via disk-based operations)

441  {
442  // Create filename from stem
443  std::string mesh_file_name =
444  Gmsh_parameters_pt->geo_and_msh_file_stem() + ".msh";
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  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
oomph::FiniteElement::size
double size() const
Definition: elements.cc:4290
oomph::Mesh::add_boundary_node
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
oomph::Mesh::Node_pt
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
oomph::Mesh::set_nboundary
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
oomph::Mesh::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
oomph::TetMeshBase::Region_attribute
Vector< double > Region_attribute
Definition: tet_mesh.h:1027
oomph::TetMeshBase::Region_element_pt
Vector< Vector< FiniteElement * > > Region_element_pt
Definition: tet_mesh.h:1022
Eigen::placeholders::end
static constexpr lastp1_t end
Definition: IndexedViewHelper.h:79
s
RealScalar s
Definition: level1_cplx_impl.h:130
calibrate.line
line
Definition: calibrate.py:103
oomph::StringConversion::to_string
std::string to_string(T object, unsigned float_precision=8)
Definition: oomph_utilities.h:189
plotDoE.x
list x
Definition: plotDoE.py:28
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
sub
EIGEN_DONT_INLINE T sub(T a, T b)
Definition: svd_common.h:238
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::Mesh::add_boundary_node(), e(), oomph::Mesh::Element_pt, Eigen::placeholders::end, oomph::Mesh::finite_element_pt(), oomph::GmshParameters::geo_and_msh_file_stem(), Gmsh_parameters_pt, i, j, calibrate::line, n, oomph::Mesh::nelement(), oomph::Mesh::Node_pt, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::TetMeshBase::Region_attribute, oomph::TetMeshBase::Region_element_pt, s, oomph::Mesh::set_nboundary(), size, oomph::FiniteElement::size(), oomph::Global_string_for_annotation::string(), sub(), oomph::StringConversion::to_string(), plotDoE::x, and oomph::Node::x().

Referenced by GmshTetScaffoldMesh().

◆ write_geo_file()

void oomph::GmshTetScaffoldMesh::write_geo_file ( const bool &  use_mesh_grading_from_file)
inlineprivate

Write geo file for gmsh.

1021  {
1022  // Transfer data from parameters
1023  TetMeshFacetedClosedSurface* outer_boundary_pt =
1024  Gmsh_parameters_pt->outer_boundary_pt();
1025 
1026  Vector<TetMeshFacetedSurface*> internal_surface_pt =
1027  Gmsh_parameters_pt->internal_surface_pt();
1028 
1029  double element_volume = Gmsh_parameters_pt->element_volume();
1030 
1031  std::string& target_size_file_name =
1032  Gmsh_parameters_pt->target_size_file_name();
1033 
1034  // Write gmsh geo file:
1035  std::string filename =
1036  Gmsh_parameters_pt->geo_and_msh_file_stem() + ".geo";
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  }
oomph::GmshParameters::internal_surface_pt
Vector< TetMeshFacetedSurface * > & internal_surface_pt()
Internal boundaries.
Definition: gmsh_tet_mesh.template.h:83
oomph::GmshParameters::element_volume
double & element_volume()
Uniform target element volume.
Definition: gmsh_tet_mesh.template.h:89
oomph::GmshParameters::target_size_file_name
std::string & target_size_file_name()
Definition: gmsh_tet_mesh.template.h:97
oomph::GmshParameters::outer_boundary_pt
TetMeshFacetedClosedSurface *& outer_boundary_pt()
Outer boundary.
Definition: gmsh_tet_mesh.template.h:77
f
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:237
int
return int(ret)+1
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
MergeRestartFiles.filename
string filename
Definition: MergeRestartFiles.py:39

References oomph::GmshParameters::element_volume(), f(), oomph::TetMeshFacet::facet_is_embedded_in_a_specified_region(), oomph::TetMeshFacetedSurface::facet_pt(), oomph::TetMeshFacetedClosedSurface::faceted_volume_represents_hole_for_gmsh(), MergeRestartFiles::filename, oomph::GmshParameters::geo_and_msh_file_stem(), Gmsh_parameters_pt, i, int(), oomph::GmshParameters::internal_surface_pt(), oomph::TetEdge::is_reversed(), j, n, oomph::TetMeshFacetedSurface::nfacet(), oomph::TetMeshFacet::nvertex(), oomph::TetMeshFacetedSurface::nvertex(), oomph::TetMeshFacet::one_based_adjacent_region_id(), oomph::TetMeshFacetedSurface::one_based_facet_boundary_id(), oomph::TetMeshFacet::one_based_region_that_facet_is_embedded_in(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::GmshParameters::outer_boundary_pt(), Eigen::bfloat16_impl::pow(), s, oomph::Global_string_for_annotation::string(), oomph::GmshParameters::target_size_file_name(), oomph::TetMeshFacet::vertex_pt(), oomph::TetMeshFacetedSurface::vertex_pt(), and oomph::TetMeshVertex::x().

Referenced by GmshTetScaffoldMesh().

Friends And Related Function Documentation

◆ GmshTetMesh

template<class ELEMENT >
friend class GmshTetMesh
friend

We're friends with the actual mesh.

Member Data Documentation

◆ Gmsh_parameters_pt

GmshParameters* oomph::GmshTetScaffoldMesh::Gmsh_parameters_pt
private

Parameters.

Referenced by create_mesh_from_msh_file(), GmshTetScaffoldMesh(), and write_geo_file().

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