Public Member Functions | Static Public Attributes | Protected Attributes | List of all members
oomph::TetMeshBase Class Reference

Base class for tet meshes (meshes made of 3D tet elements). More...

#include <tet_mesh.h>

+ Inheritance diagram for oomph::TetMeshBase:

Public Member Functions

 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...
 

Static Public Attributes

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 Attributes

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
 

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)
 
- 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)
 

Detailed Description

Base class for tet meshes (meshes made of 3D tet elements).

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////

Constructor & Destructor Documentation

◆ TetMeshBase() [1/2]

oomph::TetMeshBase::TetMeshBase ( )
inline

Constructor.

664 : Outer_boundary_pt(0) {}
oomph::TetMeshBase::Outer_boundary_pt
TetMeshFacetedClosedSurface * Outer_boundary_pt
Faceted surface that defines outer boundaries.
Definition: tet_mesh.h:1038

◆ TetMeshBase() [2/2]

oomph::TetMeshBase::TetMeshBase ( const TetMeshBase &  node)
delete

Broken copy constructor.

◆ ~TetMeshBase()

virtual oomph::TetMeshBase::~TetMeshBase ( )
inlinevirtual

Destructor (empty)

673 {}

Member Function Documentation

◆ assess_mesh_quality()

void oomph::TetMeshBase::assess_mesh_quality ( std::ofstream &  some_file)

Assess mesh quality: Ratio of max. edge length to min. height, so if it's very large it's BAAAAAD.

374  {
375  Vector<Vector<double>> edge(6);
376  for (unsigned e = 0; e < 6; e++)
377  {
378  edge[e].resize(3);
379  }
380  unsigned nel = this->nelement();
381  for (unsigned e = 0; e < nel; e++)
382  {
383  FiniteElement* fe_pt = this->finite_element_pt(e);
384  for (unsigned i = 0; i < 3; i++)
385  {
386  edge[0][i] = fe_pt->node_pt(2)->x(i) - fe_pt->node_pt(1)->x(i);
387  edge[1][i] = fe_pt->node_pt(0)->x(i) - fe_pt->node_pt(2)->x(i);
388  edge[2][i] = fe_pt->node_pt(1)->x(i) - fe_pt->node_pt(0)->x(i);
389  edge[3][i] = fe_pt->node_pt(3)->x(i) - fe_pt->node_pt(0)->x(i);
390  edge[4][i] = fe_pt->node_pt(3)->x(i) - fe_pt->node_pt(1)->x(i);
391  edge[5][i] = fe_pt->node_pt(3)->x(i) - fe_pt->node_pt(2)->x(i);
392  }
393 
394  double max_length = 0.0;
395  for (unsigned j = 0; j < 6; j++)
396  {
397  double length = 0.0;
398  for (unsigned i = 0; i < 3; i++)
399  {
400  length += edge[j][i] * edge[j][i];
401  }
402  length = sqrt(length);
403  if (length > max_length) max_length = length;
404  }
405 
406 
407  double min_height = DBL_MAX;
408  for (unsigned j = 0; j < 4; j++)
409  {
410  Vector<double> normal(3);
411  unsigned e0 = 0;
412  unsigned e1 = 0;
413  unsigned e2 = 0;
414  switch (j)
415  {
416  case 0:
417  e0 = 4;
418  e1 = 5;
419  e2 = 1;
420  break;
421 
422  case 1:
423  e0 = 1;
424  e1 = 3;
425  e2 = 2;
426  break;
427 
428  case 2:
429  e0 = 3;
430  e1 = 4;
431  e2 = 1;
432  break;
433 
434  case 3:
435  e0 = 1;
436  e1 = 2;
437  e2 = 3;
438  break;
439 
440  default:
441 
442  oomph_info << "never get here\n";
443  abort();
444  }
445 
446  normal[0] = edge[e0][1] * edge[e1][2] - edge[e0][2] * edge[e1][1];
447  normal[1] = edge[e0][2] * edge[e1][0] - edge[e0][0] * edge[e1][2];
448  normal[2] = edge[e0][0] * edge[e1][1] - edge[e0][1] * edge[e1][0];
449  double norm =
450  normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
451  double inv_norm = 1.0 / sqrt(norm);
452  normal[0] *= inv_norm;
453  normal[1] *= inv_norm;
454  normal[2] *= inv_norm;
455 
456  double height = fabs(edge[e2][0] * normal[0] + edge[e2][1] * normal[1] +
457  edge[e2][2] * normal[2]);
458 
459  if (height < min_height) min_height = height;
460  }
461 
462  double aspect_ratio = max_length / min_height;
463 
464  some_file << "ZONE N=4, E=1, F=FEPOINT, ET=TETRAHEDRON\n";
465  for (unsigned j = 0; j < 4; j++)
466  {
467  for (unsigned i = 0; i < 3; i++)
468  {
469  some_file << fe_pt->node_pt(j)->x(i) << " ";
470  }
471  some_file << aspect_ratio << std::endl;
472  }
473  some_file << "1 2 3 4" << std::endl;
474  }
475  }
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Global_Physical_Variables::height
double height(const double &x)
Height of domain.
Definition: simple_spine_channel.cc:429
WallFunction::normal
void normal(const Vector< double > &x, Vector< double > &normal)
Definition: free_surface_rotation.cc:65
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References e(), boost::multiprecision::fabs(), oomph::Mesh::finite_element_pt(), Global_Physical_Variables::height(), i, j, oomph::Mesh::nelement(), oomph::FiniteElement::node_pt(), WallFunction::normal(), oomph::oomph_info, sqrt(), and oomph::Node::x().

◆ boundary_element_in_region_pt()

FiniteElement* oomph::TetMeshBase::boundary_element_in_region_pt ( const unsigned &  b,
const unsigned &  r,
const unsigned &  e 
) const
inline

Return pointer to the e-th element adjacent to boundary b in region r.

819  {
820  // Use a constant iterator here to keep function "const" overall
821  std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
822  Boundary_region_element_pt[b].find(r);
823  if (it != Boundary_region_element_pt[b].end())
824  {
825  return (it->second)[e];
826  }
827  else
828  {
829  return 0;
830  }
831  }
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::TetMeshBase::Boundary_region_element_pt
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
Storage for elements adjacent to a boundary in a particular region.
Definition: tet_mesh.h:1031
Eigen::placeholders::end
static constexpr lastp1_t end
Definition: IndexedViewHelper.h:79
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20

References b, Boundary_region_element_pt, e(), Eigen::placeholders::end, and UniformPSDSelfTest::r.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::adapt(), and oomph::RefineableGmshTetMesh< ELEMENT >::adapt().

◆ face_index_at_boundary_in_region()

int oomph::TetMeshBase::face_index_at_boundary_in_region ( const unsigned &  b,
const unsigned &  r,
const unsigned &  e 
) const
inline

Return face index of the e-th element adjacent to boundary b in region r.

837  {
838  // Use a constant iterator here to keep function "const" overall
839  std::map<unsigned, Vector<int>>::const_iterator it =
840  Face_index_region_at_boundary[b].find(r);
841  if (it != Face_index_region_at_boundary[b].end())
842  {
843  return (it->second)[e];
844  }
845  else
846  {
847  std::ostringstream error_message;
848  error_message << "Face indices not set up for boundary " << b
849  << " in region " << r << "\n";
850  error_message << "This probably means that the boundary is not "
851  "adjacent to region\n";
852  throw OomphLibError(error_message.str(),
853  OOMPH_CURRENT_FUNCTION,
854  OOMPH_EXCEPTION_LOCATION);
855  }
856  }
oomph::TetMeshBase::Face_index_region_at_boundary
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Definition: tet_mesh.h:1035
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References b, e(), Eigen::placeholders::end, Face_index_region_at_boundary, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and UniformPSDSelfTest::r.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::adapt(), and oomph::RefineableGmshTetMesh< ELEMENT >::adapt().

◆ nboundary_element_in_region()

unsigned oomph::TetMeshBase::nboundary_element_in_region ( const unsigned &  b,
const unsigned &  r 
) const
inline

Return the number of elements adjacent to boundary b in region r.

799  {
800  // Need to use a constant iterator here to keep the function "const"
801  // Return an iterator to the appropriate entry, if we find it
802  std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
803  Boundary_region_element_pt[b].find(r);
804  if (it != Boundary_region_element_pt[b].end())
805  {
806  return (it->second).size();
807  }
808  // Otherwise there are no elements adjacent to boundary b in the region r
809  else
810  {
811  return 0;
812  }
813  }

References b, Boundary_region_element_pt, Eigen::placeholders::end, and UniformPSDSelfTest::r.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::adapt(), and oomph::RefineableGmshTetMesh< ELEMENT >::adapt().

◆ nregion()

unsigned oomph::TetMeshBase::nregion ( )
inline

Return the number of regions specified by attributes.

861  {
862  return Region_element_pt.size();
863  }
oomph::TetMeshBase::Region_element_pt
Vector< Vector< FiniteElement * > > Region_element_pt
Definition: tet_mesh.h:1022

References Region_element_pt.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::adapt(), and oomph::RefineableGmshTetMesh< ELEMENT >::adapt().

◆ nregion_element()

unsigned oomph::TetMeshBase::nregion_element ( const unsigned &  r)
inline

Return the number of elements in region r.

867  {
868  unsigned entry = 0;
869  bool found = false;
870  unsigned n = Region_attribute.size();
871  for (unsigned i = 0; i < n; i++)
872  {
873 #ifdef PARANOID
874  double diff =
875  fabs(Region_attribute[i] -
876  static_cast<double>(static_cast<unsigned>(Region_attribute[i])));
877  if (diff > 0.0)
878  {
879  std::ostringstream error_message;
880  error_message << "Region attributes should be unsigneds because we \n"
881  << "only use them to set region ids\n";
882  throw OomphLibError(error_message.str(),
883  OOMPH_CURRENT_FUNCTION,
884  OOMPH_EXCEPTION_LOCATION);
885  }
886 #endif
887  if (static_cast<unsigned>(Region_attribute[i]) == r)
888  {
889  entry = i;
890  found = true;
891  break;
892  }
893  }
894  if (found)
895  {
896  return Region_element_pt[entry].size();
897  }
898  else
899  {
900  return 0;
901  }
902  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
oomph::TetMeshBase::Region_attribute
Vector< double > Region_attribute
Definition: tet_mesh.h:1027
MergeRestartFiles.found
bool found
Definition: MergeRestartFiles.py:24

References boost::multiprecision::fabs(), MergeRestartFiles::found, i, n, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, UniformPSDSelfTest::r, Region_attribute, and Region_element_pt.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::adapt(), and oomph::RefineableGmshTetMesh< ELEMENT >::adapt().

◆ operator=()

void oomph::TetMeshBase::operator= ( const TetMeshBase &  )
delete

Broken assignment operator.

◆ region_attribute()

double oomph::TetMeshBase::region_attribute ( const unsigned &  i)
inline

Return the i-th region attribute (here only used as the (assumed to be unsigned) region id

907  {
908  return Region_attribute[i];
909  }

References i, and Region_attribute.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::adapt(), and oomph::RefineableGmshTetMesh< ELEMENT >::adapt().

◆ region_element_pt()

FiniteElement* oomph::TetMeshBase::region_element_pt ( const unsigned &  r,
const unsigned &  e 
)
inline

Return the e-th element in the r-th region.

913  {
914  unsigned entry = 0;
915  bool found = false;
916  unsigned n = Region_attribute.size();
917  for (unsigned i = 0; i < n; i++)
918  {
919 #ifdef PARANOID
920  double diff =
921  fabs(Region_attribute[i] -
922  static_cast<double>(static_cast<unsigned>(Region_attribute[i])));
923  if (diff > 0.0)
924  {
925  std::ostringstream error_message;
926  error_message << "Region attributes should be unsigneds because we \n"
927  << "only use the to set region ids\n";
928  throw OomphLibError(error_message.str(),
929  OOMPH_CURRENT_FUNCTION,
930  OOMPH_EXCEPTION_LOCATION);
931  }
932 #endif
933  if (static_cast<unsigned>(Region_attribute[i]) == r)
934  {
935  entry = i;
936  found = true;
937  break;
938  }
939  }
940  if (found)
941  {
942  return Region_element_pt[entry][e];
943  }
944  else
945  {
946  return 0;
947  }
948  }

References e(), boost::multiprecision::fabs(), MergeRestartFiles::found, i, n, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, UniformPSDSelfTest::r, Region_attribute, and Region_element_pt.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::adapt(), and oomph::RefineableGmshTetMesh< ELEMENT >::adapt().

◆ setup_boundary_coordinates() [1/4]

template<class ELEMENT >
void oomph::TetMeshBase::setup_boundary_coordinates ( const unsigned &  b)
inline

Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x-y coordinates in the plane of that boundary, with the x-axis along the line from the (lexicographically) "lower left" to the "upper right" node. The y axis is obtained by taking the cross-product of the positive x direction with the outer unit normal computed by the face elements (or its negative if switch_normal is set to true). Doc faces in output file (if it's open).

Note 1: Setup of boundary coordinates is not done if the boundary in question turns out to be nonplanar.

Note 2: If a triangular TetMeshFacet is associated with a boundary we also store the boundary coordinates of its vertices. They are needed to interpolated intrinsic coordinates of an associated GeomObject (if any) into the interior.

703  {
704  // Dummy file
705  std::ofstream some_file;
706 
707  // Don't switch the normal
708  bool switch_normal = false;
709  this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, some_file);
710  }

References b.

◆ setup_boundary_coordinates() [2/4]

template<class ELEMENT >
void oomph::TetMeshBase::setup_boundary_coordinates ( const unsigned &  b,
const bool &  switch_normal 
)
inline

Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x-y coordinates in the plane of that boundary, with the x-axis along the line from the (lexicographically) "lower left" to the "upper right" node. The y axis is obtained by taking the cross-product of the positive x direction with the outer unit normal computed by the face elements (or its negative if switch_normal is set to true). Doc faces in output file (if it's open).

Note 1: Setup of boundary coordinates is not done if the boundary in question turns out to be nonplanar.

Note 2: If a triangular TetMeshFacet is associated with a boundary we also store the boundary coordinates of its vertices. They are needed to interpolated intrinsic coordinates of an associated GeomObject (if any) into the interior. Final boolean argument allows switching of the direction of the outer unit normal.

735  {
736  // Dummy file
737  std::ofstream some_file;
738 
739  this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, some_file);
740  }

References b.

◆ setup_boundary_coordinates() [3/4]

template<class ELEMENT >
void oomph::TetMeshBase::setup_boundary_coordinates ( const unsigned &  b,
const bool &  switch_normal,
std::ofstream &  outfile 
)

Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x-y coordinates in the plane of that boundary, with the x-axis along the line from the (lexicographically) "lower left" to the "upper right" node. The y axis is obtained by taking the cross-product of the positive x direction with the outer unit normal computed by the face elements (or its negative if switch_normal is set to true). Doc faces in output file (if it's open).

Note 1: Setup of boundary coordinates is not done if the boundary in question turns out to be nonplanar.

Note 2: If a triangular TetMeshFacet is associated with a boundary we also store the boundary coordinates of its vertices. They are needed to interpolated intrinsic coordinates of an associated GeomObject (if any) into the interior. Boolean argument allows switching of the direction of the outer unit normal. Output file for doc.

/////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////// Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x-y coordinates in the plane of that boundary, with the x-axis along the line from the (lexicographically) "lower left" to the "upper right" node. The y axis is obtained by taking the cross-product of the positive x direction with the outer unit normal computed by the face elements (or its negative if switch_normal is set to true). Doc faces in output file (if it's open).

Note 1: Setup of boundary coordinates is not done if the boundary in question turns out to be nonplanar.

Note 2: If a triangular TetMeshFacet is associated with a boundary we also store the boundary coordinates of its vertices. They are needed to interpolated intrinsic coordinates of an associated GeomObject (if any) into the interior.

1095  {
1096  Node* lower_left_node_pt = 0;
1097  Node* upper_right_node_pt = 0;
1098 
1099  // Unit vector connecting lower left and upper right nodes
1100  Vector<double> b0(3);
1101 
1102  // ...and a vector normal to it
1103  Vector<double> b1(3);
1104 
1105 
1106  // Facet?
1107  TetMeshFacet* f_pt = 0;
1108  std::map<unsigned, TetMeshFacet*>::iterator it = Tet_mesh_facet_pt.find(b);
1109  if (it != Tet_mesh_facet_pt.end())
1110  {
1111  f_pt = (*it).second;
1112  }
1113 
1114  // std::cout << "Debug " << b; f_pt->output(std::cout);
1115 
1116  // Number of vertices
1117  unsigned nv = 0;
1118  if (f_pt != 0)
1119  {
1120  nv = f_pt->nvertex();
1121  }
1122 
1123  // Check for planarity and bail out if nodes are not in the same plane
1124  bool vertices_are_in_same_plane = true;
1125  for (unsigned do_for_real = 0; do_for_real < 2; do_for_real++)
1126  {
1127  // Temporary storage for face elements
1128  Vector<FiniteElement*> face_el_pt;
1129 
1130  // Loop over all elements on boundaries
1131  unsigned nel = this->nboundary_element(b);
1132 
1133  if (nel > 0)
1134  {
1135  // Loop over the bulk elements adjacent to boundary b
1136  for (unsigned e = 0; e < nel; e++)
1137  {
1138  // Get pointer to the bulk element that is adjacent to boundary b
1139  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
1140 
1141  // Find the index of the face of element e along boundary b
1142  int face_index = this->face_index_at_boundary(b, e);
1143 
1144  // Create new face element
1145  face_el_pt.push_back(
1146  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index));
1147 
1148  // Output faces?
1149  if (outfile.is_open())
1150  {
1151  face_el_pt[face_el_pt.size() - 1]->output(outfile);
1152  }
1153  }
1154 
1155  // Loop over all nodes to find the lower left and upper right ones
1156  lower_left_node_pt = this->boundary_node_pt(b, 0);
1157  upper_right_node_pt = this->boundary_node_pt(b, 0);
1158 
1159  // Loop over all nodes on the boundary
1160  unsigned nnod = this->nboundary_node(b);
1161  for (unsigned j = 0; j < nnod; j++)
1162  {
1163  // Get node
1164  Node* nod_pt = this->boundary_node_pt(b, j);
1165 
1166  // Primary criterion for lower left: Does it have a lower
1167  // z-coordinate?
1168  if (nod_pt->x(2) < lower_left_node_pt->x(2))
1169  {
1170  lower_left_node_pt = nod_pt;
1171  }
1172  // ...or is it a draw?
1173  else if (nod_pt->x(2) == lower_left_node_pt->x(2))
1174  {
1175  // If it's a draw: Does it have a lower y-coordinate?
1176  if (nod_pt->x(1) < lower_left_node_pt->x(1))
1177  {
1178  lower_left_node_pt = nod_pt;
1179  }
1180  // ...or is it a draw?
1181  else if (nod_pt->x(1) == lower_left_node_pt->x(1))
1182  {
1183  // If it's a draw: Does it have a lower x-coordinate?
1184  if (nod_pt->x(0) < lower_left_node_pt->x(0))
1185  {
1186  lower_left_node_pt = nod_pt;
1187  }
1188  }
1189  }
1190 
1191  // Primary criterion for upper right: Does it have a higher
1192  // z-coordinate?
1193  if (nod_pt->x(2) > upper_right_node_pt->x(2))
1194  {
1195  upper_right_node_pt = nod_pt;
1196  }
1197  // ...or is it a draw?
1198  else if (nod_pt->x(2) == upper_right_node_pt->x(2))
1199  {
1200  // If it's a draw: Does it have a higher y-coordinate?
1201  if (nod_pt->x(1) > upper_right_node_pt->x(1))
1202  {
1203  upper_right_node_pt = nod_pt;
1204  }
1205  // ...or is it a draw?
1206  else if (nod_pt->x(1) == upper_right_node_pt->x(1))
1207  {
1208  // If it's a draw: Does it have a higher x-coordinate?
1209  if (nod_pt->x(0) > upper_right_node_pt->x(0))
1210  {
1211  upper_right_node_pt = nod_pt;
1212  }
1213  }
1214  }
1215  }
1216 
1217  // Prepare storage for boundary coordinates
1218  Vector<double> zeta(2);
1219  /*std::cout << upper_right_node_pt->x(0) << " "
1220  << upper_right_node_pt->x(1) << " "
1221  << upper_right_node_pt->x(2) << " L ";
1222  std::cout << lower_left_node_pt->x(0) << " "
1223  << lower_left_node_pt->x(1) << " "
1224  << lower_left_node_pt->x(2) << "\n ";*/
1225 
1226 
1227  // Unit vector connecting lower left and upper right nodes
1228  b0[0] = upper_right_node_pt->x(0) - lower_left_node_pt->x(0);
1229  b0[1] = upper_right_node_pt->x(1) - lower_left_node_pt->x(1);
1230  b0[2] = upper_right_node_pt->x(2) - lower_left_node_pt->x(2);
1231 
1232  // Normalise
1233  double inv_length =
1234  1.0 / sqrt(b0[0] * b0[0] + b0[1] * b0[1] + b0[2] * b0[2]);
1235  b0[0] *= inv_length;
1236  b0[1] *= inv_length;
1237  b0[2] *= inv_length;
1238 
1239  // std::cout << "B0 ";
1240  // std::cout << b0[0] << " " << b0[1] << " " << b0[2] << "\n";
1241 
1242  // Get (outer) unit normal to first face element
1243  Vector<double> normal(3);
1244  Vector<double> s(2, 0.0);
1245  if (nv != 3)
1246  {
1247  dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[0])
1248  ->outer_unit_normal(s, normal);
1249  }
1250  else
1251  {
1252  double t1x = f_pt->vertex_pt(1)->x(0) - f_pt->vertex_pt(0)->x(0);
1253 
1254  double t1y = f_pt->vertex_pt(1)->x(1) - f_pt->vertex_pt(0)->x(1);
1255 
1256  double t1z = f_pt->vertex_pt(1)->x(2) - f_pt->vertex_pt(0)->x(2);
1257 
1258  double t2x = f_pt->vertex_pt(2)->x(0) - f_pt->vertex_pt(0)->x(0);
1259 
1260  double t2y = f_pt->vertex_pt(2)->x(1) - f_pt->vertex_pt(0)->x(1);
1261 
1262  double t2z = f_pt->vertex_pt(2)->x(2) - f_pt->vertex_pt(0)->x(2);
1263 
1264  normal[0] = t1y * t2z - t1z * t2y;
1265  normal[1] = t1z * t2x - t1x * t2z;
1266  normal[2] = t1x * t2y - t1y * t2x;
1267  double inv_length =
1268  1.0 / sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
1269  normal[2] * normal[2]);
1270  normal[0] *= inv_length;
1271  normal[1] *= inv_length;
1272  normal[2] *= inv_length;
1273  }
1274 
1275  if (switch_normal)
1276  {
1277  normal[0] = -normal[0];
1278  normal[1] = -normal[1];
1279  normal[2] = -normal[2];
1280  }
1281 
1282  // std::cout << "Normal ";
1283  // std::cout << normal[0] << " " << normal[1] << " " << normal[2] <<
1284  // "\n";
1285 
1286 
1287  // Check that all elements have the same normal
1288  for (unsigned e = 0; e < nel; e++)
1289  {
1290  // Get (outer) unit normal to face element
1291  Vector<double> my_normal(3);
1292  dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[e])
1293  ->outer_unit_normal(s, my_normal);
1294 
1295  // Dot product should be one!
1296  double dot_prod = normal[0] * my_normal[0] +
1297  normal[1] * my_normal[1] + normal[2] * my_normal[2];
1298 
1299 
1300  double error = abs(dot_prod) - 1.0;
1301  if (abs(error) > Tolerance_for_boundary_finding)
1302  {
1303  vertices_are_in_same_plane = false;
1304  }
1305  }
1306 
1307  // Bail out?
1308  if (vertices_are_in_same_plane)
1309  {
1310  // Cross-product to get second in-plane vector, normal to b0
1311  b1[0] = b0[1] * normal[2] - b0[2] * normal[1];
1312  b1[1] = b0[2] * normal[0] - b0[0] * normal[2];
1313  b1[2] = b0[0] * normal[1] - b0[1] * normal[0];
1314 
1315  // std::cout << "B1 ";
1316  // std::cout << b1[0] << " " << b1[1] << " " << b1[2] << "\n";
1317 
1318 
1319  // Assign boundary coordinates: projection onto the axes
1320  for (unsigned j = 0; j < nnod; j++)
1321  {
1322  // Get node
1323  Node* nod_pt = this->boundary_node_pt(b, j);
1324 
1325  // Difference vector to lower left corner
1326  Vector<double> delta(3);
1327  delta[0] = nod_pt->x(0) - lower_left_node_pt->x(0);
1328  delta[1] = nod_pt->x(1) - lower_left_node_pt->x(1);
1329  delta[2] = nod_pt->x(2) - lower_left_node_pt->x(2);
1330 
1331  /*std::cout << j << ": "
1332  << nod_pt->x(0) << " " << nod_pt->x(1)
1333  << " " << nod_pt->x(2) << " Delta ";
1334  std::cout << delta[0] << " " << delta[1] << " " << delta[2] << "\n";
1335  */
1336 
1337  // Project
1338  zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1339  zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1340 
1341  // Check:
1342  double error = pow(lower_left_node_pt->x(0) + zeta[0] * b0[0] +
1343  zeta[1] * b1[0] - nod_pt->x(0),
1344  2) +
1345  pow(lower_left_node_pt->x(1) + zeta[0] * b0[1] +
1346  zeta[1] * b1[1] - nod_pt->x(1),
1347  2) +
1348  pow(lower_left_node_pt->x(2) + zeta[0] * b0[2] +
1349  zeta[1] * b1[2] - nod_pt->x(2),
1350  2);
1351 
1352  if (sqrt(error) > Tolerance_for_boundary_finding)
1353  {
1354  std::ostringstream error_message;
1355  error_message
1356  << "Error in projection of boundary coordinates = "
1357  << sqrt(error) << " > Tolerance_for_boundary_finding = "
1358  << Tolerance_for_boundary_finding << std::endl
1359  << "nv = " << nv << std::endl
1360  << "Lower left: " << lower_left_node_pt->x(0) << " "
1361  << lower_left_node_pt->x(1) << " " << lower_left_node_pt->x(2)
1362  << " " << std::endl
1363  << "Upper right: " << upper_right_node_pt->x(0) << " "
1364  << upper_right_node_pt->x(1) << " " << upper_right_node_pt->x(2)
1365  << " " << std::endl
1366  << "Nodes: ";
1367  for (unsigned j = 0; j < nnod; j++)
1368  {
1369  // Get node
1370  Node* nod_pt = this->boundary_node_pt(b, j);
1371  error_message << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1372  << nod_pt->x(2) << " " << std::endl;
1373  }
1374  throw OomphLibError(error_message.str(),
1375  OOMPH_CURRENT_FUNCTION,
1376  OOMPH_EXCEPTION_LOCATION);
1377  }
1378 
1379  // Set boundary coordinate
1380  if (do_for_real == 1)
1381  {
1382  nod_pt->set_coordinates_on_boundary(b, zeta);
1383  }
1384  }
1385  } // End of if vertices are in the same plane
1386  }
1387 
1388 
1389  // Indicate that boundary coordinate has been set up
1390  if (do_for_real == 1)
1391  {
1392  Boundary_coordinate_exists[b] = true;
1393 
1394  // Facet associated with this boundary?
1395  if (f_pt != 0)
1396  {
1397  // Triangular facet: Set coordinates at vertices
1398  if (nv == 3)
1399  {
1400  Triangular_facet_vertex_boundary_coordinate[b].resize(3);
1401  for (unsigned j = 0; j < 3; j++)
1402  {
1403  // Two surface coordinates
1404  Triangular_facet_vertex_boundary_coordinate[b][j].resize(2);
1405 
1406  // Difference vector to lower left corner
1407  Vector<double> delta(3);
1408  delta[0] = f_pt->vertex_pt(j)->x(0) - lower_left_node_pt->x(0);
1409  delta[1] = f_pt->vertex_pt(j)->x(1) - lower_left_node_pt->x(1);
1410  delta[2] = f_pt->vertex_pt(j)->x(2) - lower_left_node_pt->x(2);
1411 
1412  // Project
1413  Vector<double> zeta(2);
1414  zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1415  zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1416 
1417  for (unsigned ii = 0; ii < 2; ii++)
1418  {
1419  Triangular_facet_vertex_boundary_coordinate[b][j][ii] =
1420  zeta[ii];
1421  }
1422  }
1423  }
1424  }
1425  }
1426 
1427  // Cleanup
1428  unsigned n = face_el_pt.size();
1429  for (unsigned e = 0; e < n; e++)
1430  {
1431  delete face_el_pt[e];
1432  }
1433 
1434  // Bail out?
1435  if (!vertices_are_in_same_plane)
1436  {
1437  /* oomph_info << "Vertices on boundary " << b */
1438  /* << " are not in the same plane; bailing out\n"; */
1439  return;
1440  }
1441  }
1442  }
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
oomph::Mesh::nboundary_node
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
oomph::Mesh::Boundary_coordinate_exists
std::vector< bool > Boundary_coordinate_exists
Definition: mesh.h:190
oomph::Mesh::face_index_at_boundary
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
Definition: mesh.h:896
oomph::Mesh::nboundary_element
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
oomph::Mesh::boundary_element_pt
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
oomph::Mesh::boundary_node_pt
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
oomph::TetMeshBase::Tet_mesh_facet_pt
std::map< unsigned, TetMeshFacet * > Tet_mesh_facet_pt
Definition: tet_mesh.h:1049
oomph::TetMeshBase::Tolerance_for_boundary_finding
static double Tolerance_for_boundary_finding
Definition: tet_mesh.h:677
oomph::TetMeshBase::Triangular_facet_vertex_boundary_coordinate
std::map< unsigned, Vector< Vector< double > > > Triangular_facet_vertex_boundary_coordinate
Definition: tet_mesh.h:1054
s
RealScalar s
Definition: level1_cplx_impl.h:130
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
MultiOpt.delta
int delta
Definition: MultiOpt.py:96
calibrate.error
int error
Definition: calibrate.py:297

References abs(), b, oomph::Mesh::Boundary_coordinate_exists, oomph::Mesh::boundary_element_pt(), oomph::Mesh::boundary_node_pt(), MultiOpt::delta, e(), calibrate::error, oomph::Mesh::face_index_at_boundary(), j, n, oomph::Mesh::nboundary_element(), oomph::Mesh::nboundary_node(), WallFunction::normal(), oomph::TetMeshFacet::nvertex(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Eigen::bfloat16_impl::pow(), s, oomph::Node::set_coordinates_on_boundary(), sqrt(), Tet_mesh_facet_pt, Tolerance_for_boundary_finding, Triangular_facet_vertex_boundary_coordinate, oomph::TetMeshFacet::vertex_pt(), oomph::Node::x(), oomph::TetMeshVertex::x(), and Eigen::zeta().

◆ setup_boundary_coordinates() [4/4]

template<class ELEMENT >
void oomph::TetMeshBase::setup_boundary_coordinates ( const unsigned &  b,
std::ofstream &  outfile 
)
inline

Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x-y coordinates in the plane of that boundary, with the x-axis along the line from the (lexicographically) "lower left" to the "upper right" node. The y axis is obtained by taking the cross-product of the positive x direction with the outer unit normal computed by the face elements (or its negative if switch_normal is set to true). Doc faces in output file (if it's open).

Note 1: Setup of boundary coordinates is not done if the boundary in question turns out to be nonplanar.

Note 2: If a triangular TetMeshFacet is associated with a boundary we also store the boundary coordinates of its vertices. They are needed to interpolated intrinsic coordinates of an associated GeomObject (if any) into the interior. Output file for doc.

789  {
790  // Don't switch the normal
791  bool switch_normal = false;
792  this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, outfile);
793  }

References b.

◆ setup_boundary_element_info() [1/2]

void oomph::TetMeshBase::setup_boundary_element_info ( )
inlinevirtual

Setup lookup schemes which establish which elements are located next to mesh's boundaries (wrapper to suppress doc).

Reimplemented from oomph::Mesh.

1008  {
1009  std::ofstream outfile;
1010  this->setup_boundary_element_info(outfile);
1011  }
oomph::TetMeshBase::setup_boundary_element_info
void setup_boundary_element_info()
Definition: tet_mesh.h:1007

Referenced by split_elements_in_corners().

◆ setup_boundary_element_info() [2/2]

void oomph::TetMeshBase::setup_boundary_element_info ( std::ostream &  outfile)
virtual

Setup lookup schemes which establish which elements are located next to mesh's boundaries. Doc in outfile (if it's open).

Setup lookup schemes which establish which elements are located next to which boundaries (Doc to outfile if it's open).

Reimplemented from oomph::Mesh.

101  {
102  // Should we document the output here
103  bool doc = false;
104 
105  if (outfile) doc = true;
106 
107  // Number of boundaries
108  unsigned nbound = nboundary();
109 
110  // Wipe/allocate storage for arrays
111  Boundary_element_pt.clear();
112  Face_index_at_boundary.clear();
113  Boundary_element_pt.resize(nbound);
114  Face_index_at_boundary.resize(nbound);
115  Lookup_for_elements_next_boundary_is_setup = false;
116 
117  // Temporary vector of vectors of pointers to elements on the boundaries:
118  // This is a vector to ensure UNIQUE ordering in all processors
119  Vector<Vector<FiniteElement*>> vector_of_boundary_element_pt;
120  vector_of_boundary_element_pt.resize(nbound);
121 
122  // Matrix map for working out the fixed face for elements on boundary
123  MapMatrixMixed<unsigned, FiniteElement*, int> face_identifier;
124 
125  // Loop over elements
126  //-------------------
127  unsigned nel = nelement();
128 
129 
130  // Get pointer to vector of boundaries that the
131  // node lives on
132  Vector<std::set<unsigned>*> boundaries_pt(4, 0);
133 
134  for (unsigned e = 0; e < nel; e++)
135  {
136  // Get pointer to element
137  FiniteElement* fe_pt = finite_element_pt(e);
138 
139  if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
140 
141  // Only include 3D elements! Some meshes contain interface elements too.
142  if (fe_pt->dim() == 3)
143  {
144  // Loop over the element's nodes and find out which boundaries they're
145  // on
146  // ----------------------------------------------------------------------
147  // We need only loop over the corner nodes
148  for (unsigned i = 0; i < 4; i++)
149  {
150  fe_pt->node_pt(i)->get_boundaries_pt(boundaries_pt[i]);
151  }
152 
153  // Find the common boundaries of each face
154  Vector<std::set<unsigned>> face(4);
155 
156  // NOTE: Face indices defined in Telements.h
157 
158  // Face 3 connnects points 0, 1 and 2
159  if (boundaries_pt[0] && boundaries_pt[1] && boundaries_pt[2])
160  {
161  std::set<unsigned> aux;
162 
163  std::set_intersection(
164  boundaries_pt[0]->begin(),
165  boundaries_pt[0]->end(),
166  boundaries_pt[1]->begin(),
167  boundaries_pt[1]->end(),
168  std::insert_iterator<std::set<unsigned>>(aux, aux.begin()));
169 
170  std::set_intersection(
171  aux.begin(),
172  aux.end(),
173  boundaries_pt[2]->begin(),
174  boundaries_pt[2]->end(),
175  std::insert_iterator<std::set<unsigned>>(face[3], face[3].begin()));
176  }
177 
178  // Face 2 connects points 0, 1 and 3
179  if (boundaries_pt[0] && boundaries_pt[1] && boundaries_pt[3])
180  {
181  std::set<unsigned> aux;
182 
183  std::set_intersection(
184  boundaries_pt[0]->begin(),
185  boundaries_pt[0]->end(),
186  boundaries_pt[1]->begin(),
187  boundaries_pt[1]->end(),
188  std::insert_iterator<std::set<unsigned>>(aux, aux.begin()));
189 
190  std::set_intersection(
191  aux.begin(),
192  aux.end(),
193  boundaries_pt[3]->begin(),
194  boundaries_pt[3]->end(),
195  std::insert_iterator<std::set<unsigned>>(face[2], face[2].begin()));
196  }
197 
198  // Face 1 connects points 0, 2 and 3
199  if (boundaries_pt[0] && boundaries_pt[2] && boundaries_pt[3])
200  {
201  std::set<unsigned> aux;
202 
203  std::set_intersection(
204  boundaries_pt[0]->begin(),
205  boundaries_pt[0]->end(),
206  boundaries_pt[2]->begin(),
207  boundaries_pt[2]->end(),
208  std::insert_iterator<std::set<unsigned>>(aux, aux.begin()));
209 
210  std::set_intersection(
211  aux.begin(),
212  aux.end(),
213  boundaries_pt[3]->begin(),
214  boundaries_pt[3]->end(),
215  std::insert_iterator<std::set<unsigned>>(face[1], face[1].begin()));
216  }
217 
218  // Face 0 connects points 1, 2 and 3
219  if (boundaries_pt[1] && boundaries_pt[2] && boundaries_pt[3])
220  {
221  std::set<unsigned> aux;
222 
223  std::set_intersection(
224  boundaries_pt[1]->begin(),
225  boundaries_pt[1]->end(),
226  boundaries_pt[2]->begin(),
227  boundaries_pt[2]->end(),
228  std::insert_iterator<std::set<unsigned>>(aux, aux.begin()));
229 
230  std::set_intersection(
231  aux.begin(),
232  aux.end(),
233  boundaries_pt[3]->begin(),
234  boundaries_pt[3]->end(),
235  std::insert_iterator<std::set<unsigned>>(face[0], face[0].begin()));
236  }
237 
238 
239  // We now know whether any faces lay on the boundaries
240  for (unsigned i = 0; i < 4; i++)
241  {
242  // How many boundaries are there
243  unsigned count = 0;
244 
245  // The number of the boundary
246  int boundary = -1;
247 
248  // Loop over all the members of the set and add to the count
249  // and set the boundary
250  for (std::set<unsigned>::iterator it = face[i].begin();
251  it != face[i].end();
252  ++it)
253  {
254  ++count;
255  boundary = *it;
256  }
257 
258  // If we're on more than one boundary, this is weird, so die
259  if (count > 1)
260  {
261  std::ostringstream error_stream;
262  fe_pt->output(error_stream);
263  error_stream << "Face " << i << " is on " << count
264  << " boundaries.\n";
265  error_stream << "This is rather strange.\n";
266  error_stream << "Your mesh may be too coarse or your tet mesh\n";
267  error_stream << "may be screwed up. I'm skipping the automated\n";
268  error_stream << "setup of the elements next to the boundaries\n";
269  error_stream << "lookup schemes.\n";
270  OomphLibWarning(error_stream.str(),
271  OOMPH_CURRENT_FUNCTION,
272  OOMPH_EXCEPTION_LOCATION);
273  }
274 
275  // If we have a boundary then add this to the appropriate set
276  if (boundary >= 0)
277  {
278  // Does the pointer already exits in the vector
279  Vector<FiniteElement*>::iterator b_el_it = std::find(
280  vector_of_boundary_element_pt[static_cast<unsigned>(boundary)]
281  .begin(),
282  vector_of_boundary_element_pt[static_cast<unsigned>(boundary)]
283  .end(),
284  fe_pt);
285 
286  // Only insert if we have not found it (i.e. got to the end)
287  if (b_el_it ==
288  vector_of_boundary_element_pt[static_cast<unsigned>(boundary)]
289  .end())
290  {
291  vector_of_boundary_element_pt[static_cast<unsigned>(boundary)]
292  .push_back(fe_pt);
293  }
294 
295  // Also set the fixed face
296  face_identifier(static_cast<unsigned>(boundary), fe_pt) = i;
297  }
298  }
299 
300  // Now we set the pointers to the boundary sets to zero
301  for (unsigned i = 0; i < 4; i++)
302  {
303  boundaries_pt[i] = 0;
304  }
305  }
306  }
307 
308  // Now copy everything across into permanent arrays
309  //-------------------------------------------------
310 
311  // Loop over boundaries
312  //---------------------
313  for (unsigned i = 0; i < nbound; i++)
314  {
315  // Number of elements on this boundary (currently stored in a set)
316  unsigned nel = vector_of_boundary_element_pt[i].size();
317 
318  // Allocate storage for the coordinate identifiers
319  Face_index_at_boundary[i].resize(nel);
320 
321  unsigned e_count = 0;
322  typedef Vector<FiniteElement*>::iterator IT;
323  for (IT it = vector_of_boundary_element_pt[i].begin();
324  it != vector_of_boundary_element_pt[i].end();
325  it++)
326  {
327  // Recover pointer to element
328  FiniteElement* fe_pt = *it;
329 
330  // Add to permanent storage
331  Boundary_element_pt[i].push_back(fe_pt);
332 
333  Face_index_at_boundary[i][e_count] = face_identifier(i, fe_pt);
334 
335  // Increment counter
336  e_count++;
337  }
338  }
339 
340 
341  // Doc?
342  //-----
343  if (doc)
344  {
345  // Loop over boundaries
346  for (unsigned i = 0; i < nbound; i++)
347  {
348  unsigned nel = Boundary_element_pt[i].size();
349  outfile << "Boundary: " << i << " is adjacent to " << nel << " elements"
350  << std::endl;
351 
352  // Loop over elements on given boundary
353  for (unsigned e = 0; e < nel; e++)
354  {
355  FiniteElement* fe_pt = Boundary_element_pt[i][e];
356  outfile << "Boundary element:" << fe_pt
357  << " Face index of boundary is "
358  << Face_index_at_boundary[i][e] << std::endl;
359  }
360  }
361  }
362 
363 
364  // Lookup scheme has now been setup yet
365  Lookup_for_elements_next_boundary_is_setup = true;
366  }
oomph::Mesh::Boundary_element_pt
Vector< Vector< FiniteElement * > > Boundary_element_pt
Definition: mesh.h:91
oomph::Mesh::Lookup_for_elements_next_boundary_is_setup
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:87
oomph::Mesh::Face_index_at_boundary
Vector< Vector< int > > Face_index_at_boundary
Definition: mesh.h:95
oomph::Mesh::nboundary
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827

References oomph::Mesh::Boundary_element_pt, oomph::FiniteElement::dim(), e(), Eigen::placeholders::end, oomph::Mesh::Face_index_at_boundary, oomph::Mesh::finite_element_pt(), oomph::Node::get_boundaries_pt(), i, oomph::Mesh::Lookup_for_elements_next_boundary_is_setup, oomph::Mesh::nboundary(), oomph::Mesh::nelement(), oomph::FiniteElement::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::FiniteElement::output().

◆ snap_nodes_onto_geometric_objects()

void oomph::TetMeshBase::snap_nodes_onto_geometric_objects ( )

Move the nodes on boundaries with associated GeomObjects so that they exactly coincide with the geometric object. This requires that the boundary coordinates are set up consistently.

Move the nodes on boundaries with associated Geometric Objects (if any) so that they exactly coincide with the geometric object. This requires that the boundary coordinates are set up consistently

484  {
485  // Backup in case elements get inverted
486  std::map<Node*, Vector<double>> old_nodal_posn;
487  std::map<Node*, Vector<double>> new_nodal_posn;
488  unsigned nnod = nnode();
489  for (unsigned j = 0; j < nnod; j++)
490  {
491  Node* nod_pt = node_pt(j);
492  Vector<double> x(3);
493  nod_pt->position(x);
494  old_nodal_posn[nod_pt] = x;
495  }
496 
497  // Loop over all boundaries
498  unsigned n_bound = this->nboundary();
499  for (unsigned b = 0; b < n_bound; b++)
500  {
501  bool do_it = true;
502 
503  // Accumulate reason for not snapping
504  std::stringstream reason;
505  reason << "Can't snap nodes on boundary " << b
506  << " onto geom object because: \n";
507 
508  TetMeshFacetedSurface* faceted_surface_pt = 0;
509  std::map<unsigned, TetMeshFacetedSurface*>::iterator it =
510  Tet_mesh_faceted_surface_pt.find(b);
511  if (it != Tet_mesh_faceted_surface_pt.end())
512  {
513  faceted_surface_pt = (*it).second;
514  }
515 
516  // Facet associated with this boundary?
517  if (faceted_surface_pt == 0)
518  {
519  reason << "-- no facets asssociated with boundary\n";
520  do_it = false;
521  }
522 
523  // Get geom object associated with facet
524  GeomObject* geom_obj_pt = 0;
525  if (do_it)
526  {
527  geom_obj_pt = faceted_surface_pt->geom_object_with_boundaries_pt();
528  if (geom_obj_pt == 0)
529  {
530  reason << "-- no geom object associated with boundary\n";
531  do_it = false;
532  }
533  }
534 
535  // Triangular facet?
536  if (Triangular_facet_vertex_boundary_coordinate[b].size() == 0)
537  {
538  reason << "-- facet has to be triangular and vertex coordinates have\n"
539  << " to have been set up\n";
540  do_it = false;
541  }
542 
543  // We need boundary coordinates!
544  if (!Boundary_coordinate_exists[b])
545  {
546  reason << "-- no boundary coordinates were set up\n";
547  do_it = false;
548  }
549 
550 
551  // Which facet is associated with this boundary?
552  unsigned facet_id_of_boundary = 0;
553  TetMeshFacet* f_pt = 0;
554  if (do_it)
555  {
556  unsigned nf = faceted_surface_pt->nfacet();
557  for (unsigned f = 0; f < nf; f++)
558  {
559  if ((faceted_surface_pt->one_based_facet_boundary_id(f) - 1) == b)
560  {
561  facet_id_of_boundary = f;
562  break;
563  }
564  }
565  f_pt = faceted_surface_pt->facet_pt(facet_id_of_boundary);
566 
567 
568  // Three vertices?
569  unsigned nv = f_pt->nvertex();
570  if (nv != 3)
571  {
572  reason << "-- number of facet vertices is " << nv
573  << " rather than 3\n";
574  do_it = false;
575  }
576 
577  // Have we set up zeta coordinates in geometric object?
578  if ((f_pt->vertex_pt(0)->zeta_in_geom_object().size() != 2) ||
579  (f_pt->vertex_pt(1)->zeta_in_geom_object().size() != 2) ||
580  (f_pt->vertex_pt(2)->zeta_in_geom_object().size() != 2))
581  {
582  reason << "-- no boundary coordinates were set up\n";
583  do_it = false;
584  }
585  }
586 
587 
588  // Are we ready to go?
589  if (!do_it)
590  {
591  const bool tell_us_why = false;
592  if (tell_us_why)
593  {
594  oomph_info << reason.str() << std::endl;
595  }
596  }
597  else
598  {
599  // Setup area coordinantes in triangular facet
600  double x1 = Triangular_facet_vertex_boundary_coordinate[b][0][0];
601  double y1 = Triangular_facet_vertex_boundary_coordinate[b][0][1];
602 
603  double x2 = Triangular_facet_vertex_boundary_coordinate[b][1][0];
604  double y2 = Triangular_facet_vertex_boundary_coordinate[b][1][1];
605 
606  double x3 = Triangular_facet_vertex_boundary_coordinate[b][2][0];
607  double y3 = Triangular_facet_vertex_boundary_coordinate[b][2][1];
608 
609  double detT = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3);
610 
611 
612  // Boundary coordinate (cartesian coordinates inside facet)
613  Vector<double> zeta(2);
614 
615  // Loop over all nodes on that boundary
616  const unsigned n_boundary_node = this->nboundary_node(b);
617  for (unsigned n = 0; n < n_boundary_node; ++n)
618  {
619  // Get the boundary node and coordinates
620  Node* const nod_pt = this->boundary_node_pt(b, n);
621  nod_pt->get_coordinates_on_boundary(b, zeta);
622 
623  // Now we have zeta, the cartesian boundary coordinates
624  // in the (assumed to be triangular) boundary facet; let's
625  // work out the area coordinates
626  // Notation as in
627  // https://en.wikipedia.org/wiki/Barycentric_coordinate_system
628  double s0 =
629  ((y2 - y3) * (zeta[0] - x3) + (x3 - x2) * (zeta[1] - y3)) / detT;
630  double s1 =
631  ((y3 - y1) * (zeta[0] - x3) + (x1 - x3) * (zeta[1] - y3)) / detT;
632  double s2 = 1.0 - s0 - s1;
633 
634  Vector<double> zeta_in_geom_obj(2, 0.0);
635  Vector<double> position_from_geom_obj(3, 0.0);
636 
637  // Vertex zeta coordinates
638  Vector<double> zeta_0(2);
639  zeta_0 = f_pt->vertex_pt(0)->zeta_in_geom_object();
640 
641  Vector<double> zeta_1(2);
642  zeta_1 = f_pt->vertex_pt(1)->zeta_in_geom_object();
643 
644  Vector<double> zeta_2(2);
645  zeta_2 = f_pt->vertex_pt(2)->zeta_in_geom_object();
646 
647 
648 #ifdef PARANOID
649 
650  // Compute zeta values of the vertices from parametrisation of
651  // boundaries
652  double tol = 1.0e-12;
653  Vector<double> zeta_from_boundary(2);
654  Vector<double> zeta_vertex(2);
655  for (unsigned v = 0; v < 3; v++)
656  {
657  zeta_vertex = f_pt->vertex_pt(v)->zeta_in_geom_object();
658  for (unsigned alt = 0; alt < 2; alt++)
659  {
660  switch (v)
661  {
662  case 0:
663 
664  if (alt == 0)
665  {
666  faceted_surface_pt->boundary_zeta01(
667  facet_id_of_boundary, 0.0, zeta_from_boundary);
668  }
669  else
670  {
671  faceted_surface_pt->boundary_zeta20(
672  facet_id_of_boundary, 1.0, zeta_from_boundary);
673  }
674  break;
675 
676  case 1:
677 
678  if (alt == 0)
679  {
680  faceted_surface_pt->boundary_zeta01(
681  facet_id_of_boundary, 1.0, zeta_from_boundary);
682  }
683  else
684  {
685  faceted_surface_pt->boundary_zeta12(
686  facet_id_of_boundary, 0.0, zeta_from_boundary);
687  }
688  break;
689 
690  case 2:
691 
692  if (alt == 0)
693  {
694  faceted_surface_pt->boundary_zeta12(
695  facet_id_of_boundary, 1.0, zeta_from_boundary);
696  }
697  else
698  {
699  faceted_surface_pt->boundary_zeta20(
700  facet_id_of_boundary, 0.0, zeta_from_boundary);
701  }
702  break;
703  }
704 
705  double error =
706  sqrt(pow((zeta_vertex[0] - zeta_from_boundary[0]), 2) +
707  pow((zeta_vertex[1] - zeta_from_boundary[1]), 2));
708  if (error > tol)
709  {
710  std::ostringstream error_message;
711  error_message
712  << "Error in parametrisation of boundary coordinates \n"
713  << "for vertex " << v << " [alt=" << alt << "] in facet "
714  << facet_id_of_boundary << " : \n"
715  << "zeta_vertex = [ " << zeta_vertex[0] << " "
716  << zeta_vertex[1] << " ] \n"
717  << "zeta_from_boundary = [ " << zeta_from_boundary[0]
718  << " " << zeta_from_boundary[1] << " ] \n"
719  << std::endl;
720  throw OomphLibError(error_message.str(),
721  OOMPH_CURRENT_FUNCTION,
722  OOMPH_EXCEPTION_LOCATION);
723  }
724  }
725  }
726 
727 #endif
728 
729  // Compute zeta values of the interpolation parameters
730  Vector<double> zeta_a(3, 0.0);
731  Vector<double> zeta_b(3, 0.0);
732  Vector<double> zeta_c(3, 0.0);
733  Vector<double> zeta_d(3, 0.0);
734  Vector<double> zeta_e(3, 0.0);
735  Vector<double> zeta_f(3, 0.0);
736  faceted_surface_pt->boundary_zeta01(facet_id_of_boundary, s1, zeta_a);
737  faceted_surface_pt->boundary_zeta01(
738  facet_id_of_boundary, 1.0 - s0, zeta_d);
739 
740  faceted_surface_pt->boundary_zeta12(facet_id_of_boundary, s2, zeta_c);
741  faceted_surface_pt->boundary_zeta12(
742  facet_id_of_boundary, 1.0 - s1, zeta_f);
743 
744  faceted_surface_pt->boundary_zeta20(
745  facet_id_of_boundary, 1.0 - s2, zeta_b);
746  faceted_surface_pt->boundary_zeta20(facet_id_of_boundary, s0, zeta_e);
747 
748  // Transfinite mapping
749  zeta_in_geom_obj[0] = s0 * (zeta_a[0] + zeta_b[0] - zeta_0[0]) +
750  s1 * (zeta_c[0] + zeta_d[0] - zeta_1[0]) +
751  s2 * (zeta_e[0] + zeta_f[0] - zeta_2[0]);
752  zeta_in_geom_obj[1] = s0 * (zeta_a[1] + zeta_b[1] - zeta_0[1]) +
753  s1 * (zeta_c[1] + zeta_d[1] - zeta_1[1]) +
754  s2 * (zeta_e[1] + zeta_f[1] - zeta_2[1]);
755 
756  unsigned n_tvalues =
757  1 + nod_pt->position_time_stepper_pt()->nprev_values();
758  for (unsigned t = 0; t < n_tvalues; ++t)
759  {
760  // Get the position according to the underlying geometric object
761  geom_obj_pt->position(t, zeta_in_geom_obj, position_from_geom_obj);
762 
763  // Move the node
764  for (unsigned i = 0; i < 3; i++)
765  {
766  nod_pt->x(t, i) = position_from_geom_obj[i];
767  }
768  }
769  }
770  }
771  }
772 
773  // Check if any element is inverted
774  bool some_element_is_inverted = false;
775  unsigned count = 0;
776  unsigned nel = nelement();
777  for (unsigned e = 0; e < nel; e++)
778  {
779  FiniteElement* el_pt = finite_element_pt(e);
780  bool passed = true;
781  el_pt->check_J_eulerian_at_knots(passed);
782  if (!passed)
783  {
784  some_element_is_inverted = true;
785  char filename[100];
786  std::ofstream some_file;
787  sprintf(filename, "overly_distorted_element%i.dat", count);
788  some_file.open(filename);
789  unsigned nnod_1d = el_pt->nnode_1d();
790  el_pt->output(some_file, nnod_1d);
791  some_file.close();
792 
793  // Reset to old nodal position
794  unsigned nnod = el_pt->nnode();
795  for (unsigned j = 0; j < nnod; j++)
796  {
797  Node* nod_pt = el_pt->node_pt(j);
798  Vector<double> x_current(3);
799  nod_pt->position(x_current);
800  new_nodal_posn[nod_pt] = x_current;
801  Vector<double> old_x(old_nodal_posn[nod_pt]);
802  for (unsigned i = 0; i < 3; i++)
803  {
804  nod_pt->x(i) = old_x[i];
805  }
806  }
807 
808  // Plot
809  sprintf(filename, "orig_overly_distorted_element%i.dat", count);
810  some_file.open(filename);
811  el_pt->output(some_file, nnod_1d);
812  some_file.close();
813 
814  // Reset
815  for (unsigned j = 0; j < nnod; j++)
816  {
817  Node* nod_pt = el_pt->node_pt(j);
818  for (unsigned i = 0; i < 3; i++)
819  {
820  nod_pt->x(i) = new_nodal_posn[nod_pt][i];
821  }
822  }
823 
824  // Bump
825  count++;
826  }
827  }
828  if (some_element_is_inverted)
829  {
830  std::ostringstream error_message;
831  error_message
832  << "A number of elements, namely: " << count
833  << " are inverted after snapping. Their shapes are in "
834  << " overly_distorted_element*.dat and "
835  "orig_overly_distorted_element*.dat"
836  << "Next person to get this error: Please implement a straightforward\n"
837  << "variant of one of the functors in src/mesh_smoothing to switch\n"
838  << "to harmonic mapping\n"
839  << std::endl;
840  throw OomphLibError(
841  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
842  }
843  else
844  {
845  oomph_info << "No elements are inverted after snapping. Yay!"
846  << std::endl;
847  }
848  }
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
oomph::Mesh::nnode
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
oomph::Mesh::node_pt
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
oomph::TetMeshBase::Tet_mesh_faceted_surface_pt
std::map< unsigned, TetMeshFacetedSurface * > Tet_mesh_faceted_surface_pt
Definition: tet_mesh.h:1045
f
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:237
Global_parameters::x1
Vector< double > x1(const Vector< double > &coord)
Cartesian coordinates centered at the point (0.5,1)
Definition: poisson/poisson_with_singularity/two_d_poisson.cc:86
Global_parameters::x2
Vector< double > x2(const Vector< double > &coord)
Cartesian coordinates centered at the point (1.5,1)
Definition: poisson/poisson_with_singularity/two_d_poisson.cc:102
MergeRestartFiles.filename
string filename
Definition: MergeRestartFiles.py:39
plotDoE.x
list x
Definition: plotDoE.py:28
plotPSD.t
t
Definition: plotPSD.py:36

References b, oomph::Mesh::Boundary_coordinate_exists, oomph::Mesh::boundary_node_pt(), oomph::TetMeshFacetedSurface::boundary_zeta01(), oomph::TetMeshFacetedSurface::boundary_zeta12(), oomph::TetMeshFacetedSurface::boundary_zeta20(), oomph::FiniteElement::check_J_eulerian_at_knots(), e(), calibrate::error, f(), oomph::TetMeshFacetedSurface::facet_pt(), MergeRestartFiles::filename, oomph::Mesh::finite_element_pt(), oomph::TetMeshFacetedSurface::geom_object_with_boundaries_pt(), oomph::Node::get_coordinates_on_boundary(), i, j, n, oomph::Mesh::nboundary(), oomph::Mesh::nboundary_node(), oomph::Mesh::nelement(), oomph::TetMeshFacetedSurface::nfacet(), oomph::FiniteElement::nnode(), oomph::Mesh::nnode(), oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), oomph::Mesh::node_pt(), oomph::TimeStepper::nprev_values(), oomph::TetMeshFacet::nvertex(), oomph::TetMeshFacetedSurface::one_based_facet_boundary_id(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::FiniteElement::output(), oomph::GeomObject::position(), oomph::Node::position(), oomph::Node::position_time_stepper_pt(), Eigen::bfloat16_impl::pow(), size, sqrt(), plotPSD::t, Tet_mesh_faceted_surface_pt, Triangular_facet_vertex_boundary_coordinate, v, oomph::TetMeshFacet::vertex_pt(), plotDoE::x, oomph::Node::x(), Global_parameters::x1(), Global_parameters::x2(), Eigen::zeta(), and oomph::TetMeshVertex::zeta_in_geom_object().

Referenced by oomph::GmshTetMesh< ELEMENT >::build_it(), oomph::RefineableTetgenMesh< ELEMENT >::RefineableTetgenMesh(), and oomph::TetgenMesh< ELEMENT >::TetgenMesh().

◆ snap_to_quadratic_surface() [1/2]

template<class ELEMENT >
void oomph::TetMeshBase::snap_to_quadratic_surface ( const Vector< unsigned > &  boundary_id,
const std::string &  quadratic_surface_file_name,
const bool &  switch_normal 
)
inline

Snap boundaries specified by the IDs listed in boundary_id to a quadratric surface, specified in the file quadratic_surface_file_name. This is usually used with vmtk-based meshes for which oomph-lib's xda to poly conversion code produces the files "quadratic_fsi_boundary.dat" and "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI boundary (for the fluid and the solid) and the quadratic representation of the outer boundary of the solid. When used with these files, the flag switch_normal should be set to true when calling the function for the outer boundary of the solid.

983  {
984  // Dummy doc info
985  DocInfo doc_info;
986  doc_info.disable_doc();
987  snap_to_quadratic_surface<ELEMENT>(
988  boundary_id, quadratic_surface_file_name, switch_normal, doc_info);
989  }

References oomph::DocInfo::disable_doc().

◆ snap_to_quadratic_surface() [2/2]

template<class ELEMENT >
void oomph::TetMeshBase::snap_to_quadratic_surface ( const Vector< unsigned > &  boundary_id,
const std::string &  quadratic_surface_file_name,
const bool &  switch_normal,
DocInfo &  doc_info 
)

Snap boundaries specified by the IDs listed in boundary_id to a quadratric surface, specified in the file quadratic_surface_file_name. This is usually used with vmtk-based meshes for which oomph-lib's xda to poly conversion code produces the files "quadratic_fsi_boundary.dat" and "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI boundary (for the fluid and the solid) and the quadratic representation of the outer boundary of the solid. When used with these files, the flag switch_normal should be set to true when calling the function for the outer boundary of the solid. The DocInfo object can be used to label optional output files. (Uses directory and label).

1464  {
1465  // Aux storage for processing input
1466  char dummy[101];
1467 
1468  // Prepare to doc nodes that couldn't be snapped
1469  std::ofstream the_file_non_snapped;
1470  std::string non_snapped_filename =
1471  "non_snapped_nodes_" + doc_info.label() + ".dat";
1472 
1473  // Read the number of nodes and elements (quadratic facets)
1474  std::ifstream infile(quadratic_surface_file_name.c_str(),
1475  std::ios_base::in);
1476  unsigned n_node;
1477  infile >> n_node;
1478 
1479  // Ignore rest of line
1480  infile.getline(dummy, 100);
1481 
1482  // Number of quadratic facets
1483  unsigned nel;
1484  infile >> nel;
1485 
1486  // Ignore rest of line
1487  infile.getline(dummy, 100);
1488 
1489  // Check that the number of elements matches
1490  if (nel != boundary_id.size())
1491  {
1492  std::ostringstream error_message;
1493  error_message << "Number of quadratic facets specified in "
1494  << quadratic_surface_file_name << "is: " << nel
1495  << "\nThis doesn't match the number of planar boundaries \n"
1496  << "specified in boundary_id which is "
1497  << boundary_id.size() << std::endl;
1498  throw OomphLibError(
1499  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1500  }
1501 
1502  // Temporary storage for face elements
1503  Vector<FreeStandingFaceElement<ELEMENT>*> face_el_pt;
1504 
1505  // Loop over quadratic face elements -- one for each facet
1506  for (unsigned e = 0; e < nel; e++)
1507  {
1508  face_el_pt.push_back(new FreeStandingFaceElement<ELEMENT>);
1509  }
1510 
1511 
1512  // Now build nodes
1513  unsigned n_dim = 3;
1514  unsigned n_position_type = 1;
1515  unsigned initial_n_value = 0;
1516  Vector<Node*> nod_pt(n_node);
1517  unsigned node_nmbr;
1518  std::map<unsigned, unsigned> node_number;
1519  std::ofstream node_file;
1520  for (unsigned j = 0; j < n_node; j++)
1521  {
1522  nod_pt[j] =
1523  new BoundaryNode<Node>(n_dim, n_position_type, initial_n_value);
1524  infile >> nod_pt[j]->x(0);
1525  infile >> nod_pt[j]->x(1);
1526  infile >> nod_pt[j]->x(2);
1527  infile >> node_nmbr;
1528  node_number[node_nmbr] = j;
1529  }
1530 
1531 
1532  // Now assign nodes to elements -- each element represents
1533  // distinct boundary; assign enumeration as specified by
1534  // boundary_id.
1535  for (unsigned e = 0; e < nel; e++)
1536  {
1537  unsigned nnod_el = face_el_pt[e]->nnode();
1538  unsigned j_global;
1539  for (unsigned j = 0; j < nnod_el; j++)
1540  {
1541  infile >> j_global;
1542  face_el_pt[e]->node_pt(j) = nod_pt[node_number[j_global]];
1543  face_el_pt[e]->node_pt(j)->add_to_boundary(boundary_id[e]);
1544  }
1545  face_el_pt[e]->set_boundary_number_in_bulk_mesh(boundary_id[e]);
1546  face_el_pt[e]->set_nodal_dimension(3);
1547  }
1548 
1549 
1550  // Setup boundary coordinates for each facet, using
1551  // the same strategy as for the simplex boundaries
1552  // (there's some code duplication here but it doesn't
1553  // seem worth it to rationlise this as the interface would
1554  // be pretty clunky).
1555  for (unsigned e = 0; e < nel; e++)
1556  {
1557  Vector<Vector<double>> vertex_boundary_coord(3);
1558 
1559  // Loop over all nodes to find the lower left and upper right ones
1560  Node* lower_left_node_pt = face_el_pt[e]->node_pt(0);
1561  Node* upper_right_node_pt = face_el_pt[e]->node_pt(0);
1562 
1563  // Loop over all vertex nodes
1564  for (unsigned j = 0; j < 3; j++)
1565  {
1566  // Get node
1567  Node* nod_pt = face_el_pt[e]->node_pt(j);
1568 
1569  // Primary criterion for lower left: Does it have a lower z-coordinate?
1570  if (nod_pt->x(2) < lower_left_node_pt->x(2))
1571  {
1572  lower_left_node_pt = nod_pt;
1573  }
1574  // ...or is it a draw?
1575  else if (nod_pt->x(2) == lower_left_node_pt->x(2))
1576  {
1577  // If it's a draw: Does it have a lower y-coordinate?
1578  if (nod_pt->x(1) < lower_left_node_pt->x(1))
1579  {
1580  lower_left_node_pt = nod_pt;
1581  }
1582  // ...or is it a draw?
1583  else if (nod_pt->x(1) == lower_left_node_pt->x(1))
1584  {
1585  // If it's a draw: Does it have a lower x-coordinate?
1586  if (nod_pt->x(0) < lower_left_node_pt->x(0))
1587  {
1588  lower_left_node_pt = nod_pt;
1589  }
1590  }
1591  }
1592 
1593  // Primary criterion for upper right: Does it have a higher
1594  // z-coordinate?
1595  if (nod_pt->x(2) > upper_right_node_pt->x(2))
1596  {
1597  upper_right_node_pt = nod_pt;
1598  }
1599  // ...or is it a draw?
1600  else if (nod_pt->x(2) == upper_right_node_pt->x(2))
1601  {
1602  // If it's a draw: Does it have a higher y-coordinate?
1603  if (nod_pt->x(1) > upper_right_node_pt->x(1))
1604  {
1605  upper_right_node_pt = nod_pt;
1606  }
1607  // ...or is it a draw?
1608  else if (nod_pt->x(1) == upper_right_node_pt->x(1))
1609  {
1610  // If it's a draw: Does it have a higher x-coordinate?
1611  if (nod_pt->x(0) > upper_right_node_pt->x(0))
1612  {
1613  upper_right_node_pt = nod_pt;
1614  }
1615  }
1616  }
1617  }
1618 
1619  // Prepare storage for boundary coordinates
1620  Vector<double> zeta(2);
1621 
1622  // Unit vector connecting lower left and upper right nodes
1623  Vector<double> b0(3);
1624  b0[0] = upper_right_node_pt->x(0) - lower_left_node_pt->x(0);
1625  b0[1] = upper_right_node_pt->x(1) - lower_left_node_pt->x(1);
1626  b0[2] = upper_right_node_pt->x(2) - lower_left_node_pt->x(2);
1627 
1628  // Normalise
1629  double inv_length =
1630  1.0 / sqrt(b0[0] * b0[0] + b0[1] * b0[1] + b0[2] * b0[2]);
1631  b0[0] *= inv_length;
1632  b0[1] *= inv_length;
1633  b0[2] *= inv_length;
1634 
1635  // Get (outer) unit normal to face element -- note that
1636  // with the enumeration chosen in oomph-lib's xda to poly
1637  // conversion code the sign below is correct for the outer
1638  // unit normal on the FSI interface.
1639  Vector<double> tang1(3);
1640  Vector<double> tang2(3);
1641  Vector<double> normal(3);
1642  tang1[0] =
1643  face_el_pt[e]->node_pt(1)->x(0) - face_el_pt[e]->node_pt(0)->x(0);
1644  tang1[1] =
1645  face_el_pt[e]->node_pt(1)->x(1) - face_el_pt[e]->node_pt(0)->x(1);
1646  tang1[2] =
1647  face_el_pt[e]->node_pt(1)->x(2) - face_el_pt[e]->node_pt(0)->x(2);
1648  tang2[0] =
1649  face_el_pt[e]->node_pt(2)->x(0) - face_el_pt[e]->node_pt(0)->x(0);
1650  tang2[1] =
1651  face_el_pt[e]->node_pt(2)->x(1) - face_el_pt[e]->node_pt(0)->x(1);
1652  tang2[2] =
1653  face_el_pt[e]->node_pt(2)->x(2) - face_el_pt[e]->node_pt(0)->x(2);
1654  normal[0] = -(tang1[1] * tang2[2] - tang1[2] * tang2[1]);
1655  normal[1] = -(tang1[2] * tang2[0] - tang1[0] * tang2[2]);
1656  normal[2] = -(tang1[0] * tang2[1] - tang1[1] * tang2[0]);
1657 
1658  // Normalise
1659  inv_length = 1.0 / sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
1660  normal[2] * normal[2]);
1661  normal[0] *= inv_length;
1662  normal[1] *= inv_length;
1663  normal[2] *= inv_length;
1664 
1665  // Change direction -- usually for outer boundary of solid
1666  if (switch_normal)
1667  {
1668  normal[0] = -normal[0];
1669  normal[1] = -normal[1];
1670  normal[2] = -normal[2];
1671  }
1672 
1673  // Cross-product to get second in-plane vector, normal to b0
1674  Vector<double> b1(3);
1675  b1[0] = b0[1] * normal[2] - b0[2] * normal[1];
1676  b1[1] = b0[2] * normal[0] - b0[0] * normal[2];
1677  b1[2] = b0[0] * normal[1] - b0[1] * normal[0];
1678 
1679  // Assign boundary coordinates for vertex nodes: projection onto the axes
1680  for (unsigned j = 0; j < 3; j++)
1681  {
1682  // Get node
1683  Node* nod_pt = face_el_pt[e]->node_pt(j);
1684 
1685  // Difference vector to lower left corner
1686  Vector<double> delta(3);
1687  delta[0] = nod_pt->x(0) - lower_left_node_pt->x(0);
1688  delta[1] = nod_pt->x(1) - lower_left_node_pt->x(1);
1689  delta[2] = nod_pt->x(2) - lower_left_node_pt->x(2);
1690 
1691  // Project
1692  zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1693  zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1694 
1695  vertex_boundary_coord[j].resize(2);
1696  vertex_boundary_coord[j][0] = zeta[0];
1697  vertex_boundary_coord[j][1] = zeta[1];
1698 
1699 #ifdef PARANOID
1700 
1701  // Check:
1702  double error = pow(lower_left_node_pt->x(0) + zeta[0] * b0[0] +
1703  zeta[1] * b1[0] - nod_pt->x(0),
1704  2) +
1705  pow(lower_left_node_pt->x(1) + zeta[0] * b0[1] +
1706  zeta[1] * b1[1] - nod_pt->x(1),
1707  2) +
1708  pow(lower_left_node_pt->x(2) + zeta[0] * b0[2] +
1709  zeta[1] * b1[2] - nod_pt->x(2),
1710  2);
1711 
1712  if (sqrt(error) > Tolerance_for_boundary_finding)
1713  {
1714  std::ostringstream error_message;
1715  error_message
1716  << "Error in setup of boundary coordinate " << sqrt(error)
1717  << std::endl
1718  << "exceeds tolerance specified by the static member data\n"
1719  << "TetMeshBase::Tolerance_for_boundary_finding = "
1720  << Tolerance_for_boundary_finding << std::endl
1721  << "This usually means that the boundary is not planar.\n\n"
1722  << "You can suppress this error message by recompiling \n"
1723  << "recompiling without PARANOID or by changing the tolerance.\n";
1724  throw OomphLibError(error_message.str(),
1725  OOMPH_CURRENT_FUNCTION,
1726  OOMPH_EXCEPTION_LOCATION);
1727  }
1728 #endif
1729 
1730  // Set boundary coordinate
1731  nod_pt->set_coordinates_on_boundary(boundary_id[e], zeta);
1732  }
1733 
1734  // Assign boundary coordinates of three midside nodes by linear
1735  // interpolation (corresponding to a flat facet)
1736 
1737  // Node 3 is between 0 and 1
1738  zeta[0] =
1739  0.5 * (vertex_boundary_coord[0][0] + vertex_boundary_coord[1][0]);
1740  zeta[1] =
1741  0.5 * (vertex_boundary_coord[0][1] + vertex_boundary_coord[1][1]);
1742  face_el_pt[e]->node_pt(3)->set_coordinates_on_boundary(boundary_id[e],
1743  zeta);
1744 
1745  // Node 4 is between 1 and 2
1746  zeta[0] =
1747  0.5 * (vertex_boundary_coord[1][0] + vertex_boundary_coord[2][0]);
1748  zeta[1] =
1749  0.5 * (vertex_boundary_coord[1][1] + vertex_boundary_coord[2][1]);
1750  face_el_pt[e]->node_pt(4)->set_coordinates_on_boundary(boundary_id[e],
1751  zeta);
1752 
1753  // Node 5 is between 2 and 0
1754  zeta[0] =
1755  0.5 * (vertex_boundary_coord[2][0] + vertex_boundary_coord[0][0]);
1756  zeta[1] =
1757  0.5 * (vertex_boundary_coord[2][1] + vertex_boundary_coord[0][1]);
1758  face_el_pt[e]->node_pt(5)->set_coordinates_on_boundary(boundary_id[e],
1759  zeta);
1760  }
1761 
1762 
1763  // Loop over elements/facets = boundaries to snap
1764  bool success = true;
1765  for (unsigned b = 0; b < nel; b++)
1766  {
1767  // Doc boundary coordinates on quadratic patches
1768  std::ofstream the_file;
1769  std::ofstream the_file_before;
1770  std::ofstream the_file_after;
1771  if (doc_info.is_doc_enabled())
1772  {
1773  std::ostringstream filename;
1774  filename << doc_info.directory() << "/quadratic_coordinates_"
1775  << doc_info.label() << b << ".dat";
1776  the_file.open(filename.str().c_str());
1777 
1778  std::ostringstream filename1;
1779  filename1 << doc_info.directory() << "/quadratic_nodes_before_"
1780  << doc_info.label() << b << ".dat";
1781  the_file_before.open(filename1.str().c_str());
1782 
1783  std::ostringstream filename2;
1784  filename2 << doc_info.directory() << "/quadratic_nodes_after_"
1785  << doc_info.label() << b << ".dat";
1786  the_file_after.open(filename2.str().c_str());
1787  }
1788 
1789  // Cast the element pointer
1790  FreeStandingFaceElement<ELEMENT>* el_pt = face_el_pt[b];
1791 
1792  // Doc boundary coordinate on quadratic facet representation
1793  if (doc_info.is_doc_enabled())
1794  {
1795  Vector<double> s(2);
1796  Vector<double> zeta(2);
1797  Vector<double> x(3);
1798  unsigned n_plot = 3;
1799  the_file << el_pt->tecplot_zone_string(n_plot);
1800 
1801  // Loop over plot points
1802  unsigned num_plot_points = el_pt->nplot_points(n_plot);
1803  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1804  {
1805  // Get local coordinates of plot point
1806  el_pt->get_s_plot(iplot, n_plot, s);
1807  el_pt->interpolated_zeta(s, zeta);
1808  el_pt->interpolated_x(s, x);
1809  for (unsigned i = 0; i < 3; i++)
1810  {
1811  the_file << x[i] << " ";
1812  }
1813  for (unsigned i = 0; i < 2; i++)
1814  {
1815  the_file << zeta[i] << " ";
1816  }
1817  the_file << std::endl;
1818  }
1819  el_pt->write_tecplot_zone_footer(the_file, n_plot);
1820 
1821  // std::cout << "Facet " << b << " corresponds to quadratic
1822  // boundary "
1823  // << boundary_id[b] << " which contains "
1824  // << this->nboundary_element(boundary_id[b])
1825  // << " element[s] " << std::endl;
1826  }
1827 
1828 
1829  // Loop over bulk elements that are adjacent to quadratic boundary
1830  Vector<double> boundary_zeta(2);
1831  Vector<double> quadratic_coordinates(3);
1832  GeomObject* geom_obj_pt = 0;
1833  Vector<double> s_geom_obj(2);
1834  unsigned nel = this->nboundary_element(boundary_id[b]);
1835  for (unsigned e = 0; e < nel; e++)
1836  {
1837  // Get pointer to the bulk element that is adjacent to boundary
1838  FiniteElement* bulk_elem_pt =
1839  this->boundary_element_pt(boundary_id[b], e);
1840 
1841  // Loop over nodes
1842  unsigned nnod = bulk_elem_pt->nnode();
1843  for (unsigned j = 0; j < nnod; j++)
1844  {
1845  Node* nod_pt = bulk_elem_pt->node_pt(j);
1846  if (nod_pt->is_on_boundary(boundary_id[b]))
1847  {
1848  nod_pt->get_coordinates_on_boundary(boundary_id[b], boundary_zeta);
1849 
1850  // Doc it?
1851  if (doc_info.is_doc_enabled())
1852  {
1853  the_file_before << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1854  << nod_pt->x(2) << " " << boundary_zeta[0] << " "
1855  << boundary_zeta[1] << " " << std::endl;
1856  }
1857 
1858  // Find local coordinate in quadratic facet
1859  el_pt->locate_zeta(boundary_zeta, geom_obj_pt, s_geom_obj);
1860  if (geom_obj_pt != 0)
1861  {
1862  geom_obj_pt->position(s_geom_obj, quadratic_coordinates);
1863  nod_pt->x(0) = quadratic_coordinates[0];
1864  nod_pt->x(1) = quadratic_coordinates[1];
1865  nod_pt->x(2) = quadratic_coordinates[2];
1866  }
1867  else
1868  {
1869  // Get ready to bail out below...
1870  success = false;
1871 
1872  std::ostringstream error_message;
1873  error_message
1874  << "Warning: Couldn't find GeomObject during snapping to\n"
1875  << "quadratic surface for boundary " << boundary_id[b]
1876  << ". I'm leaving the node where it was. Will bail out "
1877  "later.\n";
1878  OomphLibWarning(error_message.str(),
1879  "TetgenMesh::snap_to_quadratic_surface()",
1880  OOMPH_EXCEPTION_LOCATION);
1881  if (!the_file_non_snapped.is_open())
1882  {
1883  the_file_non_snapped.open(non_snapped_filename.c_str());
1884  }
1885  the_file_non_snapped << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1886  << nod_pt->x(2) << " " << boundary_zeta[0]
1887  << " " << boundary_zeta[1] << " "
1888  << std::endl;
1889  }
1890 
1891  // Doc it?
1892  if (doc_info.is_doc_enabled())
1893  {
1894  the_file_after << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1895  << nod_pt->x(2) << " " << boundary_zeta[0] << " "
1896  << boundary_zeta[1] << " " << std::endl;
1897  }
1898  }
1899  }
1900  }
1901 
1902  // Close doc file
1903  the_file.close();
1904  the_file_before.close();
1905  the_file_after.close();
1906  }
1907 
1908  // Bail out?
1909  if (!success)
1910  {
1911  std::ostringstream error_message;
1912  error_message
1913  << "Warning: Couldn't find GeomObject during snapping to\n"
1914  << "quadratic surface. Bailing out.\n"
1915  << "Nodes that couldn't be snapped are contained in \n"
1916  << "file: " << non_snapped_filename << ".\n"
1917  << "This problem may arise because the switch_normal flag was \n"
1918  << "set wrongly.\n";
1919  throw OomphLibError(
1920  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1921  }
1922 
1923  // Close
1924  if (!the_file_non_snapped.is_open())
1925  {
1926  the_file_non_snapped.close();
1927  }
1928 
1929  // Kill auxiliary FreeStandingFaceElements
1930  for (unsigned e = 0; e < nel; e++)
1931  {
1932  delete face_el_pt[e];
1933  face_el_pt[e] = 0;
1934  }
1935 
1936  // Kill boundary nodes
1937  unsigned nn = nod_pt.size();
1938  for (unsigned j = 0; j < nn; j++)
1939  {
1940  delete nod_pt[j];
1941  }
1942  }
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

References b, oomph::Mesh::boundary_element_pt(), MultiOpt::delta, oomph::DocInfo::directory(), e(), calibrate::error, MergeRestartFiles::filename, oomph::Node::get_coordinates_on_boundary(), i, oomph::DocInfo::is_doc_enabled(), oomph::Node::is_on_boundary(), j, oomph::DocInfo::label(), oomph::Mesh::nboundary_element(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), WallFunction::normal(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::GeomObject::position(), Eigen::bfloat16_impl::pow(), s, oomph::Node::set_coordinates_on_boundary(), sqrt(), oomph::Global_string_for_annotation::string(), Tolerance_for_boundary_finding, plotDoE::x, oomph::Node::x(), and Eigen::zeta().

◆ split_elements_in_corners()

template<class ELEMENT >
void oomph::TetMeshBase::split_elements_in_corners ( TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper)

Non-Delaunay split elements that have three faces on a boundary into sons. Timestepper species timestepper for new nodes; defaults to to steady timestepper.

Non-delaunay split elements that have three faces on a boundary into sons.

1951  {
1952  // Setup boundary lookup scheme if required
1953  if (!Lookup_for_elements_next_boundary_is_setup)
1954  {
1955  setup_boundary_element_info();
1956  }
1957 
1958  // Find out how many nodes we have along each element edge
1959  unsigned nnode_1d = finite_element_pt(0)->nnode_1d();
1960  // Find out the total number of nodes
1961  unsigned nnode = this->finite_element_pt(0)->nnode();
1962 
1963  // At the moment we're only able to deal with nnode_1d=2 or 3.
1964  if ((nnode_1d != 2) && (nnode_1d != 3))
1965  {
1966  std::ostringstream error_message;
1967  error_message << "Mesh generation from tetgen currently only works\n";
1968  error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
1969  error_message << "for nnode_1d=" << nnode_1d << std::endl;
1970 
1971  throw OomphLibError(
1972  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1973  }
1974 
1975  // Map to store how many boundaries elements are located on
1976  std::map<FiniteElement*, unsigned> count;
1977 
1978  // Loop over all boundaries
1979  unsigned nb = this->nboundary();
1980  for (unsigned b = 0; b < nb; b++)
1981  {
1982  // Loop over elements next to that boundary
1983  unsigned nel = this->nboundary_element(b);
1984  for (unsigned e = 0; e < nel; e++)
1985  {
1986  // Get pointer to element
1987  FiniteElement* el_pt = boundary_element_pt(b, e);
1988 
1989  // Bump up counter
1990  count[el_pt]++;
1991  }
1992  }
1993 
1994  // To avoid having to check the map for all elements (which will
1995  // inflate it to the size of all elements!), move offending elements
1996  // into set
1997  std::set<FiniteElement*> elements_to_be_split;
1998  for (std::map<FiniteElement*, unsigned>::iterator it = count.begin();
1999  it != count.end();
2000  it++)
2001  {
2002  // Get pointer to element: Key
2003  FiniteElement* el_pt = it->first;
2004 
2005  // Does it have to be split?
2006  if (it->second > 2)
2007  {
2008  elements_to_be_split.insert(el_pt);
2009  }
2010  }
2011 
2012  // Vector for retained or newly built elements
2013  unsigned nel = this->nelement();
2014  Vector<FiniteElement*> new_or_retained_el_pt;
2015  new_or_retained_el_pt.reserve(nel);
2016 
2017  // Map which returns the 4 newly created elements for each old corner
2018  // element
2019  std::map<FiniteElement*, Vector<FiniteElement*>>
2020  old_to_new_corner_element_map;
2021 
2022  // Now loop over all elements
2023  for (unsigned e = 0; e < nel; e++)
2024  {
2025  // Get pointer to element
2026  FiniteElement* el_pt = this->finite_element_pt(e);
2027 
2028  // Does it have to be split? I.e. is it in the set?
2029  std::set<FiniteElement*>::iterator it = std::find(
2030  elements_to_be_split.begin(), elements_to_be_split.end(), el_pt);
2031 
2032  // It's not in the set, so iterator has reached end
2033  if (it == elements_to_be_split.end())
2034  {
2035  // Carry it across
2036  new_or_retained_el_pt.push_back(el_pt);
2037  }
2038  // It's in the set of elements to be split
2039  else
2040  {
2041  // Storage for new nodes -- Note: All new nodes are interior and
2042  // therefore cannot be boundary nodes!
2043  Node* node0_pt = 0;
2044  Node* node1_pt = 0;
2045  Node* node2_pt = 0;
2046  Node* node3_pt = 0;
2047  Node* node4_pt = 0;
2048  Node* node5_pt = 0;
2049  Node* node6_pt = 0;
2050  Node* node7_pt = 0;
2051  Node* node8_pt = 0;
2052  Node* node9_pt = 0;
2053  Node* node10_pt = 0;
2054 
2055  // Create first new element
2056  FiniteElement* el1_pt = new ELEMENT;
2057 
2058  // Copy existing nodes
2059  el1_pt->node_pt(0) = el_pt->node_pt(0);
2060  el1_pt->node_pt(1) = el_pt->node_pt(1);
2061  el1_pt->node_pt(3) = el_pt->node_pt(3);
2062  if (nnode_1d == 3)
2063  {
2064  el1_pt->node_pt(4) = el_pt->node_pt(4);
2065  el1_pt->node_pt(6) = el_pt->node_pt(6);
2066  el1_pt->node_pt(9) = el_pt->node_pt(9);
2067  }
2068 
2069  // Create new nodes
2070  // If we have an enriched element then don't need to construct centroid
2071  // node
2072  if (nnode == 15)
2073  {
2074  node0_pt = el_pt->node_pt(14);
2075  el1_pt->node_pt(2) = node0_pt;
2076  node5_pt = el1_pt->construct_node(11, time_stepper_pt);
2077  node6_pt = el1_pt->construct_node(13, time_stepper_pt);
2078  node9_pt = el1_pt->construct_node(12, time_stepper_pt);
2079 
2080  // Copy others over
2081  el1_pt->node_pt(10) = el_pt->node_pt(10);
2082  }
2083  // If not enriched we do
2084  else
2085  {
2086  node0_pt = el1_pt->construct_node(2, time_stepper_pt);
2087  }
2088  if (nnode_1d == 3)
2089  {
2090  node1_pt = el1_pt->construct_boundary_node(5, time_stepper_pt);
2091  node2_pt = el1_pt->construct_boundary_node(7, time_stepper_pt);
2092  node4_pt = el1_pt->construct_boundary_node(8, time_stepper_pt);
2093  }
2094 
2095 
2096  // Create second new element
2097  FiniteElement* el2_pt = new ELEMENT;
2098 
2099  // Copy existing nodes
2100  el2_pt->node_pt(0) = el_pt->node_pt(0);
2101  el2_pt->node_pt(1) = el_pt->node_pt(1);
2102  el2_pt->node_pt(2) = el_pt->node_pt(2);
2103  if (nnode_1d == 3)
2104  {
2105  el2_pt->node_pt(4) = el_pt->node_pt(4);
2106  el2_pt->node_pt(5) = el_pt->node_pt(5);
2107  el2_pt->node_pt(7) = el_pt->node_pt(7);
2108  }
2109 
2110  // Create new node
2111  if (nnode_1d == 3)
2112  {
2113  node3_pt = el2_pt->construct_boundary_node(8, time_stepper_pt);
2114  }
2115 
2116  // Copy existing new nodes
2117  el2_pt->node_pt(3) = node0_pt;
2118  if (nnode_1d == 3)
2119  {
2120  el2_pt->node_pt(6) = node1_pt;
2121  el2_pt->node_pt(9) = node2_pt;
2122  }
2123 
2124  // Copy and constuct other nodes for enriched elements
2125  if (nnode == 15)
2126  {
2127  el2_pt->node_pt(10) = node5_pt;
2128  el2_pt->node_pt(11) = el_pt->node_pt(11);
2129  node8_pt = el2_pt->construct_node(12, time_stepper_pt);
2130  node10_pt = el2_pt->construct_node(13, time_stepper_pt);
2131  }
2132 
2133  // Create third new element
2134  FiniteElement* el3_pt = new ELEMENT;
2135 
2136  // Copy existing nodes
2137  el3_pt->node_pt(1) = el_pt->node_pt(1);
2138  el3_pt->node_pt(2) = el_pt->node_pt(2);
2139  el3_pt->node_pt(3) = el_pt->node_pt(3);
2140  if (nnode_1d == 3)
2141  {
2142  el3_pt->node_pt(7) = el_pt->node_pt(7);
2143  el3_pt->node_pt(8) = el_pt->node_pt(8);
2144  el3_pt->node_pt(9) = el_pt->node_pt(9);
2145  }
2146 
2147  // Copy existing new nodes
2148  el3_pt->node_pt(0) = node0_pt;
2149  if (nnode_1d == 3)
2150  {
2151  el3_pt->node_pt(4) = node2_pt;
2152  el3_pt->node_pt(5) = node3_pt;
2153  el3_pt->node_pt(6) = node4_pt;
2154  }
2155 
2156  // Copy and constuct other nodes for enriched elements
2157  if (nnode == 15)
2158  {
2159  el3_pt->node_pt(10) = node6_pt;
2160  el3_pt->node_pt(11) = node10_pt;
2161  node7_pt = el3_pt->construct_node(12, time_stepper_pt);
2162  el3_pt->node_pt(13) = el_pt->node_pt(13);
2163  }
2164 
2165 
2166  // Create fourth new element
2167  FiniteElement* el4_pt = new ELEMENT;
2168 
2169  // Copy existing nodes
2170  el4_pt->node_pt(0) = el_pt->node_pt(0);
2171  el4_pt->node_pt(2) = el_pt->node_pt(2);
2172  el4_pt->node_pt(3) = el_pt->node_pt(3);
2173  if (nnode_1d == 3)
2174  {
2175  el4_pt->node_pt(5) = el_pt->node_pt(5);
2176  el4_pt->node_pt(6) = el_pt->node_pt(6);
2177  el4_pt->node_pt(8) = el_pt->node_pt(8);
2178  }
2179 
2180  // Copy existing new nodes
2181  el4_pt->node_pt(1) = node0_pt;
2182  if (nnode_1d == 3)
2183  {
2184  el4_pt->node_pt(4) = node1_pt;
2185  el4_pt->node_pt(7) = node3_pt;
2186  el4_pt->node_pt(9) = node4_pt;
2187  }
2188 
2189  // Copy all other nodes
2190  if (nnode == 15)
2191  {
2192  el4_pt->node_pt(10) = node9_pt;
2193  el4_pt->node_pt(11) = node8_pt;
2194  el4_pt->node_pt(12) = el_pt->node_pt(12);
2195  el4_pt->node_pt(13) = node7_pt;
2196  ;
2197  }
2198 
2199 
2200  // Add new elements and nodes
2201  new_or_retained_el_pt.push_back(el1_pt);
2202  new_or_retained_el_pt.push_back(el2_pt);
2203  new_or_retained_el_pt.push_back(el3_pt);
2204  new_or_retained_el_pt.push_back(el4_pt);
2205 
2206  // create a vector of the newly created elements
2207  Vector<FiniteElement*> temp_vec(4);
2208  temp_vec[0] = el1_pt;
2209  temp_vec[1] = el2_pt;
2210  temp_vec[2] = el3_pt;
2211  temp_vec[3] = el4_pt;
2212 
2213  // add the vector to the map
2214  old_to_new_corner_element_map.insert(
2215  std::pair<FiniteElement*, Vector<FiniteElement*>>(el_pt, temp_vec));
2216 
2217  if (nnode != 15)
2218  {
2219  this->add_node_pt(node0_pt);
2220  }
2221  this->add_node_pt(node1_pt);
2222  this->add_node_pt(node2_pt);
2223  this->add_node_pt(node3_pt);
2224  this->add_node_pt(node4_pt);
2225  if (nnode == 15)
2226  {
2227  this->add_node_pt(node5_pt);
2228  this->add_node_pt(node6_pt);
2229  this->add_node_pt(node7_pt);
2230  this->add_node_pt(node8_pt);
2231  this->add_node_pt(node9_pt);
2232  }
2233 
2234  // Set nodal positions
2235  for (unsigned i = 0; i < 3; i++)
2236  {
2237  // Only bother to set centroid if does not already exist
2238  if (nnode != 15)
2239  {
2240  node0_pt->x(i) =
2241  0.25 * (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i) +
2242  el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i));
2243  }
2244 
2245  if (nnode_1d == 3)
2246  {
2247  node1_pt->x(i) = 0.5 * (el_pt->node_pt(0)->x(i) + node0_pt->x(i));
2248  node2_pt->x(i) = 0.5 * (el_pt->node_pt(1)->x(i) + node0_pt->x(i));
2249  node3_pt->x(i) = 0.5 * (el_pt->node_pt(2)->x(i) + node0_pt->x(i));
2250  node4_pt->x(i) = 0.5 * (el_pt->node_pt(3)->x(i) + node0_pt->x(i));
2251  }
2252  }
2253 
2254 
2255  // Construct the four interior nodes if needed
2256  // and add to the list
2257  if (nnode == 15)
2258  {
2259  // Set the positions of the newly created mid-face nodes
2260  // New Node 5 lies in the plane between original nodes 0 1 centroid
2261  for (unsigned i = 0; i < 3; ++i)
2262  {
2263  node5_pt->x(i) =
2264  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i) +
2265  el_pt->node_pt(14)->x(i)) /
2266  3.0;
2267  }
2268 
2269  // New Node 6 lies in the plane between original nodes 1 3 centroid
2270  for (unsigned i = 0; i < 3; ++i)
2271  {
2272  node6_pt->x(i) =
2273  (el_pt->node_pt(1)->x(i) + el_pt->node_pt(3)->x(i) +
2274  el_pt->node_pt(14)->x(i)) /
2275  3.0;
2276  }
2277 
2278  // New Node 7 lies in the plane between original nodes 2 3 centroid
2279  for (unsigned i = 0; i < 3; ++i)
2280  {
2281  node7_pt->x(i) =
2282  (el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i) +
2283  el_pt->node_pt(14)->x(i)) /
2284  3.0;
2285  }
2286 
2287  // New Node 8 lies in the plane between original nodes 0 2 centroid
2288  for (unsigned i = 0; i < 3; ++i)
2289  {
2290  node8_pt->x(i) =
2291  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(2)->x(i) +
2292  el_pt->node_pt(14)->x(i)) /
2293  3.0;
2294  }
2295 
2296  // New Node 9 lies in the plane between original nodes 0 3 centroid
2297  for (unsigned i = 0; i < 3; ++i)
2298  {
2299  node9_pt->x(i) =
2300  (el_pt->node_pt(0)->x(i) + el_pt->node_pt(3)->x(i) +
2301  el_pt->node_pt(14)->x(i)) /
2302  3.0;
2303  }
2304 
2305  // New Node 10 lies in the plane between original nodes 1 2 centroid
2306  for (unsigned i = 0; i < 3; ++i)
2307  {
2308  node10_pt->x(i) =
2309  (el_pt->node_pt(1)->x(i) + el_pt->node_pt(2)->x(i) +
2310  el_pt->node_pt(14)->x(i)) /
2311  3.0;
2312  }
2313 
2314  // Now create the new centroid nodes
2315 
2316  // First element
2317  Node* temp_nod_pt = el1_pt->construct_node(14, time_stepper_pt);
2318  for (unsigned i = 0; i < 3; ++i)
2319  {
2320  double av_pos = 0.0;
2321  for (unsigned j = 0; j < 4; j++)
2322  {
2323  av_pos += el1_pt->node_pt(j)->x(i);
2324  }
2325 
2326  temp_nod_pt->x(i) = 0.25 * av_pos;
2327  }
2328 
2329  this->add_node_pt(temp_nod_pt);
2330 
2331  // Second element
2332  temp_nod_pt = el2_pt->construct_node(14, time_stepper_pt);
2333  for (unsigned i = 0; i < 3; ++i)
2334  {
2335  double av_pos = 0.0;
2336  for (unsigned j = 0; j < 4; j++)
2337  {
2338  av_pos += el2_pt->node_pt(j)->x(i);
2339  }
2340  temp_nod_pt->x(i) = 0.25 * av_pos;
2341  }
2342  this->add_node_pt(temp_nod_pt);
2343 
2344  // Third element
2345  temp_nod_pt = el3_pt->construct_node(14, time_stepper_pt);
2346  for (unsigned i = 0; i < 3; ++i)
2347  {
2348  double av_pos = 0.0;
2349  for (unsigned j = 0; j < 4; j++)
2350  {
2351  av_pos += el3_pt->node_pt(j)->x(i);
2352  }
2353  temp_nod_pt->x(i) = 0.25 * av_pos;
2354  }
2355  this->add_node_pt(temp_nod_pt);
2356 
2357  // Fourth element
2358  temp_nod_pt = el4_pt->construct_node(14, time_stepper_pt);
2359  for (unsigned i = 0; i < 3; ++i)
2360  {
2361  double av_pos = 0.0;
2362  for (unsigned j = 0; j < 4; j++)
2363  {
2364  av_pos += el4_pt->node_pt(j)->x(i);
2365  }
2366  temp_nod_pt->x(i) = 0.25 * av_pos;
2367  }
2368  this->add_node_pt(temp_nod_pt);
2369  }
2370 
2371  // Kill the old element
2372  delete el_pt;
2373  }
2374  }
2375 
2376  // Flush element storage
2377  Element_pt.clear();
2378 
2379  // Copy across
2380  nel = new_or_retained_el_pt.size();
2381  Element_pt.resize(nel);
2382  for (unsigned e = 0; e < nel; e++)
2383  {
2384  Element_pt[e] = new_or_retained_el_pt[e];
2385  }
2386 
2387  // Setup boundary lookup scheme again
2388  setup_boundary_element_info();
2389 
2390  // -------------------------------------------------------------------------
2391  // The various boundary/region lookups now need updating to account for the
2392  // newly added/removed elements. This will be done in two stages:
2393  // Step 1: Add the new elements to the vector of elements in the same region
2394  // as the original corner element, and then delete the originals.
2395  // Updating this lookup makes things easier in the following step.
2396  // Step 2: Regenerate the two more specific lookups: One which gives the
2397  // elements on a given boundary in a given region, and the other
2398  // which maps elements on a given boundary in a given region to the
2399  // element's face index on that boundary.
2400  //
2401  // N.B. the lookup Triangular_facet_vertex_boundary_coordinate is setup in
2402  // the call to setup_boundary_element_info() above so doesn't need
2403  // additional work.
2404 
2405  // if we have no regions then we have no lookups to update so we're done
2406  // here
2407  if (Region_attribute.size() == 0)
2408  {
2409  return;
2410  }
2411  // if we haven't had to split any corner elements then don't need to fiddle
2412  // with the lookups
2413  if (old_to_new_corner_element_map.size() == 0)
2414  {
2415  oomph_info << "\nNo corner elements need splitting\n\n";
2416  return;
2417  }
2418 
2419  // ------------------------------------------
2420  // Step 1: update the region element lookup
2421 
2422  // loop over the map of old corner elements which have been split
2423  for (std::map<FiniteElement*, Vector<FiniteElement*>>::iterator map_it =
2424  old_to_new_corner_element_map.begin();
2425  map_it != old_to_new_corner_element_map.end();
2426  map_it++)
2427  {
2428  // extract the old and new elements from the map
2429  FiniteElement* original_el_pt = map_it->first;
2430  Vector<FiniteElement*> new_el_pt = map_it->second;
2431 
2432  unsigned original_region_index = 0;
2433 
2434 #ifdef PARANOID
2435  // flag for paranoia, if for some reason we don't find the original corner
2436  // element in any of the regions
2437  bool found = false;
2438 #endif
2439 
2440  Vector<FiniteElement*>::iterator region_element_it;
2441 
2442  // loop over the regions and look for this original corner element to find
2443  // out which region it used to be in, so we can add the new elements to
2444  // the same region.
2445  for (unsigned region_index = 0; region_index < Region_element_pt.size();
2446  region_index++)
2447  {
2448  // for each region, search the vector of elements in this region for the
2449  // original corner element
2450  region_element_it = std::find(Region_element_pt[region_index].begin(),
2451  Region_element_pt[region_index].end(),
2452  original_el_pt);
2453 
2454  // if the iterator hasn't reached the end then we've found it
2455  if (region_element_it != Region_element_pt[region_index].end())
2456  {
2457  // save the region index we're at
2458  original_region_index = region_index;
2459 
2460 #ifdef PARANOID
2461  // update the paranoid flag
2462  found = true;
2463 #endif
2464 
2465  // regions can't overlap, so once we've found one we're done
2466  break;
2467  }
2468  }
2469 
2470 #ifdef PARANOID
2471  if (!found)
2472  {
2473  std::ostringstream error_message;
2474  error_message
2475  << "The corner element being split does not appear to be in any "
2476  << "region, so something has gone wrong with the region lookup "
2477  "scheme\n";
2478 
2479  throw OomphLibError(error_message.str(),
2480  OOMPH_CURRENT_FUNCTION,
2481  OOMPH_EXCEPTION_LOCATION);
2482  }
2483 #endif
2484 
2485  // Now update the basic region lookup. The iterator will still point to
2486  // the original corner element in the lookup, so we can delete this easily
2487  Region_element_pt[original_region_index].erase(region_element_it);
2488  for (unsigned i = 0; i < 4; i++)
2489  {
2490  Region_element_pt[original_region_index].push_back(new_el_pt[i]);
2491  }
2492  }
2493  // ------------------------------------------
2494  // Step 2: Clear and regenerate lookups
2495 
2496  Face_index_region_at_boundary.clear();
2497  Boundary_region_element_pt.clear();
2498 
2499  Face_index_region_at_boundary.resize(nboundary());
2500  Boundary_region_element_pt.resize(nboundary());
2501 
2502  for (unsigned b = 0; b < nboundary(); b++)
2503  {
2504  // Loop over elements next to that boundary
2505  unsigned nel = this->nboundary_element(b);
2506  for (unsigned e = 0; e < nel; e++)
2507  {
2508  FiniteElement* el_pt = boundary_element_pt(b, e);
2509 
2510  // now search for it in each region
2511  for (unsigned r_index = 0; r_index < Region_attribute.size(); r_index++)
2512  {
2513  unsigned region_id = static_cast<unsigned>(Region_attribute[r_index]);
2514 
2515  Vector<FiniteElement*>::iterator it =
2516  std::find(Region_element_pt[r_index].begin(),
2517  Region_element_pt[r_index].end(),
2518  el_pt);
2519 
2520  // if we find this element in the current region, then update our
2521  // lookups
2522  if (it != Region_element_pt[r_index].end())
2523  {
2524  Boundary_region_element_pt[b][region_id].push_back(el_pt);
2525 
2526  unsigned face_index = face_index_at_boundary(b, e);
2527  Face_index_region_at_boundary[b][region_id].push_back(face_index);
2528  }
2529  }
2530  }
2531  }
2532 
2533  oomph_info << "\nNumber of outer corner elements split: "
2534  << old_to_new_corner_element_map.size() << "\n\n";
2535 
2536  } // end split_elements_in_corners()
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::FiniteElement::nnode_1d
virtual unsigned nnode_1d() const
Definition: elements.h:2218
oomph::Mesh::add_node_pt
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
oomph::Mesh::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
nb
int nb
Definition: level2_impl.h:286

References oomph::Mesh::add_node_pt(), b, oomph::Mesh::boundary_element_pt(), Boundary_region_element_pt, oomph::FiniteElement::construct_boundary_node(), oomph::FiniteElement::construct_node(), e(), oomph::Mesh::Element_pt, Eigen::placeholders::end, oomph::Mesh::face_index_at_boundary(), Face_index_region_at_boundary, oomph::Mesh::finite_element_pt(), MergeRestartFiles::found, i, j, oomph::Mesh::Lookup_for_elements_next_boundary_is_setup, nb, oomph::Mesh::nboundary(), oomph::Mesh::nboundary_element(), oomph::Mesh::nelement(), oomph::FiniteElement::nnode(), oomph::Mesh::nnode(), oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Region_attribute, Region_element_pt, setup_boundary_element_info(), and oomph::Node::x().

Member Data Documentation

◆ Boundary_region_element_pt

Vector<std::map<unsigned, Vector<FiniteElement*> > > oomph::TetMeshBase::Boundary_region_element_pt
protected

Storage for elements adjacent to a boundary in a particular region.

Referenced by boundary_element_in_region_pt(), oomph::GmshTetMesh< ELEMENT >::build_from_scaffold(), nboundary_element_in_region(), and split_elements_in_corners().

◆ Face_index_region_at_boundary

Vector<std::map<unsigned, Vector<int> > > oomph::TetMeshBase::Face_index_region_at_boundary
protected

Storage for the face index adjacent to a boundary in a particular region

Referenced by oomph::GmshTetMesh< ELEMENT >::build_from_scaffold(), face_index_at_boundary_in_region(), and split_elements_in_corners().

◆ Internal_surface_pt

Vector<TetMeshFacetedSurface*> oomph::TetMeshBase::Internal_surface_pt
protected

Vector to faceted surfaces that define internal boundaries.

Referenced by oomph::GmshTetMesh< ELEMENT >::build_it(), oomph::RefineableTetgenMesh< ELEMENT >::RefineableTetgenMesh(), and oomph::TetgenMesh< ELEMENT >::TetgenMesh().

◆ Outer_boundary_pt

TetMeshFacetedClosedSurface* oomph::TetMeshBase::Outer_boundary_pt
protected

Faceted surface that defines outer boundaries.

Referenced by oomph::GmshTetMesh< ELEMENT >::build_it(), oomph::RefineableTetgenMesh< ELEMENT >::RefineableTetgenMesh(), and oomph::TetgenMesh< ELEMENT >::TetgenMesh().

◆ Region_attribute

Vector<double> oomph::TetMeshBase::Region_attribute
protected

Vector of attributes associated with the elements in each region NOTE: double is enforced on us by tetgen. We use it as an unsigned to indicate the actual (zero-based) region ID

Referenced by oomph::GmshTetMesh< ELEMENT >::build_from_scaffold(), oomph::GmshTetScaffoldMesh::create_mesh_from_msh_file(), nregion_element(), region_attribute(), region_element_pt(), and split_elements_in_corners().

◆ Region_element_pt

Vector<Vector<FiniteElement*> > oomph::TetMeshBase::Region_element_pt
protected

Vectors of vectors of elements in each region (note: this just stores them; the region IDs are contained in Region_attribute!)

Referenced by oomph::GmshTetMesh< ELEMENT >::build_from_scaffold(), oomph::GmshTetScaffoldMesh::create_mesh_from_msh_file(), nregion(), nregion_element(), region_element_pt(), and split_elements_in_corners().

◆ Tet_mesh_facet_pt

std::map<unsigned, TetMeshFacet*> oomph::TetMeshBase::Tet_mesh_facet_pt
protected

Reverse lookup scheme: Pointer to facet (if any!) associated with boundary b

Referenced by oomph::GmshTetMesh< ELEMENT >::build_it(), and setup_boundary_coordinates().

◆ Tet_mesh_faceted_surface_pt

std::map<unsigned, TetMeshFacetedSurface*> oomph::TetMeshBase::Tet_mesh_faceted_surface_pt
protected

Reverse lookup scheme: Pointer to faceted surface (if any!) associated with boundary b

Referenced by oomph::GmshTetMesh< ELEMENT >::build_it(), and snap_nodes_onto_geometric_objects().

◆ Time_stepper_pt

TimeStepper* oomph::TetMeshBase::Time_stepper_pt
protected

Timestepper used to build nodes.

Referenced by oomph::GmshTetMesh< ELEMENT >::build_it(), oomph::RefineableTetgenMesh< ELEMENT >::RefineableTetgenMesh(), oomph::TetgenMesh< ELEMENT >::set_mesh_level_time_stepper(), and oomph::TetgenMesh< ELEMENT >::TetgenMesh().

◆ Tolerance_for_boundary_finding

double oomph::TetMeshBase::Tolerance_for_boundary_finding = 1.0e-5
static

Global static data that specifies the permitted error in the setup of the boundary coordinates

Referenced by setup_boundary_coordinates(), and snap_to_quadratic_surface().

◆ Triangular_facet_vertex_boundary_coordinate

std::map<unsigned, Vector<Vector<double> > > oomph::TetMeshBase::Triangular_facet_vertex_boundary_coordinate
protected

Boundary coordinates of vertices in triangular facets associated with given boundary. Is only set up for triangular facets!

Referenced by setup_boundary_coordinates(), and snap_nodes_onto_geometric_objects().

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