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

#include <unstructured_two_d_mesh_geometry_base.h>

+ Inheritance diagram for oomph::UnstructuredTwoDMeshGeometryBase:

Public Member Functions

 UnstructuredTwoDMeshGeometryBase ()
 Empty constructor. More...
 
 UnstructuredTwoDMeshGeometryBase (const UnstructuredTwoDMeshGeometryBase &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const UnstructuredTwoDMeshGeometryBase &)=delete
 Broken assignment operator. More...
 
 ~UnstructuredTwoDMeshGeometryBase ()
 Empty destructor. More...
 
unsigned nregion ()
 Return the number of regions specified by attributes. More...
 
unsigned nregion_element (const unsigned &i)
 Return the number of elements in the i-th region. More...
 
FiniteElement * region_element_pt (const unsigned &i, const unsigned &e)
 Return the e-th element in the i-th region. More...
 
unsigned nregion_attribute ()
 Return the number of attributes used in the mesh. More...
 
double region_attribute (const unsigned &i)
 Return the attribute associated with region i. More...
 
GeomObject * boundary_geom_object_pt (const unsigned &b)
 
std::map< unsigned, GeomObject * > & boundary_geom_object_pt ()
 Return direct access to the geometric object storage. More...
 
std::map< unsigned, Vector< double > > & boundary_coordinate_limits ()
 
Vector< double > & boundary_coordinate_limits (const unsigned &b)
 
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...
 
TriangleMeshPolyLine * boundary_polyline_pt (const unsigned &b)
 
std::map< unsigned, std::set< Node * > > & nodes_on_boundary_pt ()
 
const bool get_connected_vertex_number_on_destination_polyline (TriangleMeshPolyLine *dst_polyline_pt, Vector< double > &vertex_coordinates, unsigned &vertex_number)
 
void check_contiguousness_on_polylines_helper (Vector< TriangleMeshPolyLine * > &polylines_pt, unsigned &index)
 
void check_contiguousness_on_polylines_helper (Vector< TriangleMeshPolyLine * > &polylines_pt, unsigned &index_halo_start, unsigned &index_halo_end)
 
bool is_point_inside_polygon_helper (Vector< Vector< double >> polygon_vertices, Vector< double > point)
 Helper function that checks if a given point is inside a polygon. More...
 
void enable_automatic_creation_of_vertices_on_boundaries ()
 
void disable_automatic_creation_of_vertices_on_boundaries ()
 
bool is_automatic_creation_of_vertices_on_boundaries_allowed ()
 
template<class ELEMENT >
void setup_boundary_coordinates (const unsigned &b)
 
template<class ELEMENT >
void setup_boundary_coordinates (const unsigned &b, std::ofstream &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 setup_boundary_element_info ()
 
virtual void setup_boundary_element_info (std::ostream &outfile)
 
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 bool Suppress_warning_about_regions_and_boundaries
 Public static flag to suppress warning; defaults to false. More...
 
- Static Public Attributes inherited from oomph::Mesh
static Steady< 0 > Default_TimeStepper
 The Steady Timestepper. More...
 
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
 Static boolean flag to control warning about mesh level timesteppers. More...
 

Protected Member Functions

void snap_nodes_onto_geometric_objects ()
 
void copy_connection_information (TriangleMeshCurveSection *input_curve_pt, TriangleMeshCurveSection *output_curve_pt)
 
void copy_connection_information_to_sub_polylines (TriangleMeshCurveSection *input_curve_pt, TriangleMeshCurveSection *output_curve_pt)
 
- Protected Member Functions inherited from oomph::Mesh
unsigned long assign_global_eqn_numbers (Vector< double * > &Dof_pt)
 Assign (global) equation numbers to the nodes. More...
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign local equation numbers in all elements. More...
 
void convert_to_boundary_node (Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
 
void convert_to_boundary_node (Node *&node_pt)
 

Protected Attributes

bool Allow_automatic_creation_of_vertices_on_boundaries
 
std::map< unsigned, Vector< FiniteElement * > > Region_element_pt
 
Vector< double > Region_attribute
 Vector of attributes associated with the elements in each region. More...
 
std::map< unsigned, GeomObject * > Boundary_geom_object_pt
 Storage for the geometric objects associated with any boundaries. More...
 
std::map< unsigned, Vector< double > > Boundary_coordinate_limits
 
Vector< TriangleMeshPolygon * > Outer_boundary_pt
 Polygon that defines outer boundaries. More...
 
Vector< TriangleMeshPolygon * > Internal_polygon_pt
 Vector of polygons that define internal polygons. More...
 
Vector< TriangleMeshOpenCurve * > Internal_open_curve_pt
 Vector of open polylines that define internal curves. More...
 
Vector< Vector< double > > Extra_holes_coordinates
 Storage for extra coordinates for holes. More...
 
std::map< unsigned, Vector< double > > Regions_coordinates
 
std::map< unsigned, TriangleMeshCurveSection * > Boundary_curve_section_pt
 
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
 Storage for the face index adjacent to a boundary in a particular region. More...
 
std::map< unsigned, Vector< std::pair< double, double > > > Polygonal_vertex_arclength_info
 
std::map< unsigned, std::set< Node * > > Nodes_on_boundary_pt
 
std::set< TriangleMeshCurveSection * > Free_curve_section_pt
 
std::set< TriangleMeshPolygon * > Free_polygon_pt
 
std::set< TriangleMeshOpenCurve * > Free_open_curve_pt
 
- 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)
 

Detailed Description

Contains functions which define the geometry of the mesh, i.e. regions, boundaries, etc.

Constructor & Destructor Documentation

◆ UnstructuredTwoDMeshGeometryBase() [1/2]

oomph::UnstructuredTwoDMeshGeometryBase::UnstructuredTwoDMeshGeometryBase ( )
inline

Empty constructor.

1744 {}

◆ UnstructuredTwoDMeshGeometryBase() [2/2]

oomph::UnstructuredTwoDMeshGeometryBase::UnstructuredTwoDMeshGeometryBase ( const UnstructuredTwoDMeshGeometryBase &  dummy)
delete

Broken copy constructor.

◆ ~UnstructuredTwoDMeshGeometryBase()

oomph::UnstructuredTwoDMeshGeometryBase::~UnstructuredTwoDMeshGeometryBase ( )
inline

Empty destructor.

1754 {}

Member Function Documentation

◆ boundary_coordinate_limits() [1/2]

std::map<unsigned, Vector<double> >& oomph::UnstructuredTwoDMeshGeometryBase::boundary_coordinate_limits ( )
inline

Return access to the vector of boundary coordinates associated with each geometric object

1850  {
1851  return Boundary_coordinate_limits;
1852  }
oomph::UnstructuredTwoDMeshGeometryBase::Boundary_coordinate_limits
std::map< unsigned, Vector< double > > Boundary_coordinate_limits
Definition: unstructured_two_d_mesh_geometry_base.h:2593

References Boundary_coordinate_limits.

Referenced by setup_boundary_coordinates().

◆ boundary_coordinate_limits() [2/2]

Vector<double>& oomph::UnstructuredTwoDMeshGeometryBase::boundary_coordinate_limits ( const unsigned &  b)
inline

Return access to the coordinate limits associated with the geometric object associated with boundary b

1857  {
1858  std::map<unsigned, Vector<double>>::iterator it;
1859  it = Boundary_coordinate_limits.find(b);
1860  if (it == Boundary_coordinate_limits.end())
1861  {
1862  throw OomphLibError(
1863  "No coordinate limits associated with this boundary\n",
1864  OOMPH_CURRENT_FUNCTION,
1865  OOMPH_EXCEPTION_LOCATION);
1866  }
1867  else
1868  {
1869  return (*it).second;
1870  }
1871  }
b
Scalar * b
Definition: benchVecAdd.cpp:17
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, Boundary_coordinate_limits, OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ boundary_element_in_region_pt()

FiniteElement* oomph::UnstructuredTwoDMeshGeometryBase::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.

1896  {
1897  // Use a constant iterator here to keep function "const" overall
1898  std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
1899  Boundary_region_element_pt[b].find(r);
1900  if (it != Boundary_region_element_pt[b].end())
1901  {
1902  return (it->second)[e];
1903  }
1904  else
1905  {
1906  return 0;
1907  }
1908  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::UnstructuredTwoDMeshGeometryBase::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: unstructured_two_d_mesh_geometry_base.h:2617
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 setup_boundary_coordinates().

◆ boundary_geom_object_pt() [1/2]

std::map<unsigned, GeomObject*>& oomph::UnstructuredTwoDMeshGeometryBase::boundary_geom_object_pt ( )
inline

Return direct access to the geometric object storage.

1843  {
1844  return Boundary_geom_object_pt;
1845  }
oomph::UnstructuredTwoDMeshGeometryBase::Boundary_geom_object_pt
std::map< unsigned, GeomObject * > Boundary_geom_object_pt
Storage for the geometric objects associated with any boundaries.
Definition: unstructured_two_d_mesh_geometry_base.h:2589

References Boundary_geom_object_pt.

Referenced by setup_boundary_coordinates(), and snap_nodes_onto_geometric_objects().

◆ boundary_geom_object_pt() [2/2]

GeomObject* oomph::UnstructuredTwoDMeshGeometryBase::boundary_geom_object_pt ( const unsigned &  b)
inline

Return the geometric object associated with the b-th boundary or null if the boundary has associated geometric object.

1828  {
1829  std::map<unsigned, GeomObject*>::iterator it;
1830  it = Boundary_geom_object_pt.find(b);
1831  if (it == Boundary_geom_object_pt.end())
1832  {
1833  return 0;
1834  }
1835  else
1836  {
1837  return (*it).second;
1838  }
1839  }

References b, and Boundary_geom_object_pt.

◆ boundary_polyline_pt()

TriangleMeshPolyLine* oomph::UnstructuredTwoDMeshGeometryBase::boundary_polyline_pt ( const unsigned &  b)
inline

Return pointer to the current polyline that describes the b-th mesh boundary

1938  {
1939  std::map<unsigned, TriangleMeshCurveSection*>::iterator it =
1940  Boundary_curve_section_pt.find(b);
1941  // Search for the polyline associated with the given boundary
1942  if (it != Boundary_curve_section_pt.end())
1943  {
1944  return dynamic_cast<TriangleMeshPolyLine*>(
1945  Boundary_curve_section_pt[b]);
1946  }
1947  // If the boundary was not established then return 0, or a null pointer
1948  return 0;
1949  }
oomph::UnstructuredTwoDMeshGeometryBase::Boundary_curve_section_pt
std::map< unsigned, TriangleMeshCurveSection * > Boundary_curve_section_pt
Definition: unstructured_two_d_mesh_geometry_base.h:2613

References b, and Boundary_curve_section_pt.

Referenced by ContactProblem< ELEMENT >::actions_before_adapt().

◆ check_contiguousness_on_polylines_helper() [1/2]

void oomph::UnstructuredTwoDMeshGeometryBase::check_contiguousness_on_polylines_helper ( Vector< TriangleMeshPolyLine * > &  polylines_pt,
unsigned &  index 
)

Sort the polylines coming from joining them. Check whether it is necessary to reverse them or not. Used when joining two curve polylines in order to create a polygon

◆ check_contiguousness_on_polylines_helper() [2/2]

void oomph::UnstructuredTwoDMeshGeometryBase::check_contiguousness_on_polylines_helper ( Vector< TriangleMeshPolyLine * > &  polylines_pt,
unsigned &  index_halo_start,
unsigned &  index_halo_end 
)

Sort the polylines coming from joining them. Check whether it is necessary to reverse them or not. Used when joining polylines and they still do not create a polygon

◆ copy_connection_information()

void oomph::UnstructuredTwoDMeshGeometryBase::copy_connection_information ( TriangleMeshCurveSection *  input_curve_pt,
TriangleMeshCurveSection *  output_curve_pt 
)
protected

Helper function to copy the connection information from the input curve(polyline or curviline) to the output polyline

◆ copy_connection_information_to_sub_polylines()

void oomph::UnstructuredTwoDMeshGeometryBase::copy_connection_information_to_sub_polylines ( TriangleMeshCurveSection *  input_curve_pt,
TriangleMeshCurveSection *  output_curve_pt 
)
protected

Helper function to copy the connection information from the input curve(polyline or curviline) to the output sub-polyline

◆ disable_automatic_creation_of_vertices_on_boundaries()

void oomph::UnstructuredTwoDMeshGeometryBase::disable_automatic_creation_of_vertices_on_boundaries ( )
inline

Disables the creation of points (by Triangle) on the outer and internal boundaries

1994  {
1995  Allow_automatic_creation_of_vertices_on_boundaries = false;
1996  }
oomph::UnstructuredTwoDMeshGeometryBase::Allow_automatic_creation_of_vertices_on_boundaries
bool Allow_automatic_creation_of_vertices_on_boundaries
Definition: unstructured_two_d_mesh_geometry_base.h:2575

References Allow_automatic_creation_of_vertices_on_boundaries.

◆ enable_automatic_creation_of_vertices_on_boundaries()

void oomph::UnstructuredTwoDMeshGeometryBase::enable_automatic_creation_of_vertices_on_boundaries ( )
inline

Enables the creation of points (by Triangle) on the outer and internal boundaries

1987  {
1988  Allow_automatic_creation_of_vertices_on_boundaries = true;
1989  }

References Allow_automatic_creation_of_vertices_on_boundaries.

◆ face_index_at_boundary_in_region()

int oomph::UnstructuredTwoDMeshGeometryBase::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.

1914  {
1915  // Use a constant iterator here to keep function "const" overall
1916  std::map<unsigned, Vector<int>>::const_iterator it =
1917  Face_index_region_at_boundary[b].find(r);
1918  if (it != Face_index_region_at_boundary[b].end())
1919  {
1920  return (it->second)[e];
1921  }
1922  else
1923  {
1924  std::ostringstream error_message;
1925  error_message << "Face indices not set up for boundary " << b
1926  << " in region " << r << "\n";
1927  error_message << "This probably means that the boundary is not "
1928  "adjacent to region\n";
1929  throw OomphLibError(error_message.str(),
1930  OOMPH_CURRENT_FUNCTION,
1931  OOMPH_EXCEPTION_LOCATION);
1932  }
1933  }
oomph::UnstructuredTwoDMeshGeometryBase::Face_index_region_at_boundary
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Storage for the face index adjacent to a boundary in a particular region.
Definition: unstructured_two_d_mesh_geometry_base.h:2620

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

Referenced by setup_boundary_coordinates().

◆ get_connected_vertex_number_on_destination_polyline()

const bool oomph::UnstructuredTwoDMeshGeometryBase::get_connected_vertex_number_on_destination_polyline ( TriangleMeshPolyLine *  dst_polyline_pt,
Vector< double > &  vertex_coordinates,
unsigned &  vertex_number 
)

Gets the vertex number on the destination polyline (used to create the connections among shared boundaries)

Gets the vertex number on the destination polyline (used / to create the connections among shared boundaries)

4125  {
4126  // The returning flag indicating that the vertex was found in the
4127  // destination boundary
4128  bool found_vertex_number = false;
4129 
4130  // Get the number of vertices in the destination polyline
4131  const unsigned n_vertices = dst_polyline_pt->nvertex();
4132 
4133  // Loop over the vertices and return the closest vertex
4134  // to the given vertex coordinates
4135  for (unsigned i = 0; i < n_vertices; i++)
4136  {
4137  // Get the i-vertex on the polyline
4138  Vector<double> current_vertex = dst_polyline_pt->vertex_coordinate(i);
4139 
4140  // Compute the distance to the input vertex
4141  double dist = (vertex_coordinates[0] - current_vertex[0]) *
4142  (vertex_coordinates[0] - current_vertex[0]) +
4143  (vertex_coordinates[1] - current_vertex[1]) *
4144  (vertex_coordinates[1] - current_vertex[1]);
4145  dist = sqrt(dist);
4146 
4147  // Have we found the vertex?
4148  if (dist < ToleranceForVertexMismatchInPolygons::Tolerable_error)
4149  {
4150  // Set the vertex index
4151  vertex_number = i;
4152 
4153  // Set the flag to found
4154  found_vertex_number = true;
4155 
4156  // Break the serching
4157  break;
4158  }
4159  } // for (i < n_vertices)
4160 
4161  // Return if the connection was found
4162  return found_vertex_number;
4163  }
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
oomph::ToleranceForVertexMismatchInPolygons::Tolerable_error
double Tolerable_error
Definition: unstructured_two_d_mesh_geometry_base.cc:1103

References i, oomph::TriangleMeshPolyLine::nvertex(), sqrt(), oomph::ToleranceForVertexMismatchInPolygons::Tolerable_error, and oomph::TriangleMeshPolyLine::vertex_coordinate().

◆ is_automatic_creation_of_vertices_on_boundaries_allowed()

bool oomph::UnstructuredTwoDMeshGeometryBase::is_automatic_creation_of_vertices_on_boundaries_allowed ( )
inline

Returns the status of the variable Allow_automatic_creation_of_vertices_on_boundaries

2001  {
2002  return Allow_automatic_creation_of_vertices_on_boundaries;
2003  }

References Allow_automatic_creation_of_vertices_on_boundaries.

◆ is_point_inside_polygon_helper()

bool oomph::UnstructuredTwoDMeshGeometryBase::is_point_inside_polygon_helper ( Vector< Vector< double >>  polygon_vertices,
Vector< double >  point 
)

Helper function that checks if a given point is inside a polygon.

Helper function that checks if a given point is inside a polygon (a set of sorted vertices that connected create a polygon)

4067  {
4068  // Total number of vertices (the first and last vertex should be the same)
4069  const unsigned n_vertex = polygon_vertices.size();
4070 
4071  // Check if internal point is actually located in bounding polygon
4072  // Reference: http://paulbourke.net/geometry/insidepoly/
4073 
4074  // Counter for number of intersections
4075  unsigned intersect_counter = 0;
4076 
4077  // Get first vertex
4078  Vector<double> p1 = polygon_vertices[0];
4079  for (unsigned i = 1; i <= n_vertex; i++)
4080  {
4081  // Get second vertex by wrap-around
4082  Vector<double> p2 = polygon_vertices[i % n_vertex];
4083 
4084  if (point[1] > std::min(p1[1], p2[1]))
4085  {
4086  if (point[1] <= std::max(p1[1], p2[1]))
4087  {
4088  if (point[0] <= std::max(p1[0], p2[0]))
4089  {
4090  if (p1[1] != p2[1])
4091  {
4092  double xintersect =
4093  (point[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0];
4094  if ((p1[0] == p2[0]) || (point[0] <= xintersect))
4095  {
4096  intersect_counter++;
4097  }
4098  }
4099  }
4100  }
4101  }
4102  p1 = p2;
4103  }
4104 
4105  // Even number of intersections: outside
4106  if (intersect_counter % 2 == 0)
4107  {
4108  return false;
4109  }
4110  else
4111  {
4112  return true;
4113  }
4114  }
p1
Vector3f p1
Definition: MatrixBase_all.cpp:2
min
#define min(a, b)
Definition: datatypes.h:22
max
#define max(a, b)
Definition: datatypes.h:23

References i, max, min, and p1.

◆ nboundary_element_in_region()

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

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

1876  {
1877  // Need to use a constant iterator here to keep the function "const".
1878  // Return an iterator to the appropriate entry, if we find it
1879  std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
1880  Boundary_region_element_pt[b].find(r);
1881  if (it != Boundary_region_element_pt[b].end())
1882  {
1883  return (it->second).size();
1884  }
1885  // Otherwise there are no elements adjacent to boundary b in the region r
1886  else
1887  {
1888  return 0;
1889  }
1890  }

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

Referenced by setup_boundary_coordinates().

◆ nodes_on_boundary_pt()

std::map<unsigned, std::set<Node*> >& oomph::UnstructuredTwoDMeshGeometryBase::nodes_on_boundary_pt ( )
inline

Gets a pointer to a set with all the nodes related with a boundary

1954  {
1955  return Nodes_on_boundary_pt;
1956  }
oomph::UnstructuredTwoDMeshGeometryBase::Nodes_on_boundary_pt
std::map< unsigned, std::set< Node * > > Nodes_on_boundary_pt
Definition: unstructured_two_d_mesh_geometry_base.h:2634

References Nodes_on_boundary_pt.

◆ nregion()

unsigned oomph::UnstructuredTwoDMeshGeometryBase::nregion ( )
inline

Return the number of regions specified by attributes.

1758  {
1759  return Region_element_pt.size();
1760  }
oomph::UnstructuredTwoDMeshGeometryBase::Region_element_pt
std::map< unsigned, Vector< FiniteElement * > > Region_element_pt
Definition: unstructured_two_d_mesh_geometry_base.h:2583

References Region_element_pt.

Referenced by setup_boundary_coordinates().

◆ nregion_attribute()

unsigned oomph::UnstructuredTwoDMeshGeometryBase::nregion_attribute ( )
inline

Return the number of attributes used in the mesh.

1815  {
1816  return Region_attribute.size();
1817  }
oomph::UnstructuredTwoDMeshGeometryBase::Region_attribute
Vector< double > Region_attribute
Vector of attributes associated with the elements in each region.
Definition: unstructured_two_d_mesh_geometry_base.h:2586

References Region_attribute.

◆ nregion_element()

unsigned oomph::UnstructuredTwoDMeshGeometryBase::nregion_element ( const unsigned &  i)
inline

Return the number of elements in the i-th region.

1764  {
1765  // Create an iterator to iterate over Region_element_pt
1766  std::map<unsigned, Vector<FiniteElement*>>::iterator it;
1767 
1768  // Find the entry of Region_element_pt associated with the i-th region
1769  it = Region_element_pt.find(i);
1770 
1771  // If there is an entry associated with the i-th region
1772  if (it != Region_element_pt.end())
1773  {
1774  return (it->second).size();
1775  }
1776  else
1777  {
1778  return 0;
1779  }
1780  } // End of nregion_element

References i, and Region_element_pt.

◆ operator=()

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

Broken assignment operator.

◆ region_attribute()

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

Return the attribute associated with region i.

1821  {
1822  return Region_attribute[i];
1823  }

References i, and Region_attribute.

◆ region_element_pt()

FiniteElement* oomph::UnstructuredTwoDMeshGeometryBase::region_element_pt ( const unsigned &  i,
const unsigned &  e 
)
inline

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

1784  {
1785  // Create an iterator to iterate over Region_element_pt
1786  std::map<unsigned, Vector<FiniteElement*>>::iterator it;
1787 
1788  // Find the entry of Region_element_pt associated with the i-th region
1789  it = Region_element_pt.find(i);
1790 
1791  // If there is an entry associated with the i-th region
1792  if (it != Region_element_pt.end())
1793  {
1794  // Return a pointer to the e-th element in the i-th region
1795  return (it->second)[e];
1796  }
1797  else
1798  {
1799  // Create a stringstream object
1800  std::stringstream error_message;
1801 
1802  // Create the error message
1803  error_message << "There are no regions associated with "
1804  << "region ID " << i << ".";
1805 
1806  // Throw an error
1807  throw OomphLibError(error_message.str(),
1808  OOMPH_CURRENT_FUNCTION,
1809  OOMPH_EXCEPTION_LOCATION);
1810  }
1811  } // End of region_element_pt

References e(), i, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Region_element_pt.

◆ setup_boundary_coordinates() [1/2]

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

Setup boundary coordinate on boundary b. Boundary coordinate increases continously along polygonal boundary. It's zero at the lowest left node on the boundary.

2443  {
2444  // Dummy file
2445  std::ofstream some_file;
2446  setup_boundary_coordinates<ELEMENT>(b, some_file);
2447  }

References b.

◆ setup_boundary_coordinates() [2/2]

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

Setup boundary coordinate on boundary b. Doc Faces in outfile. Boundary coordinate increases continously along polygonal boundary. It's zero at the lowest left node on the boundary.

Setup boundary coordinate on boundary b. Doc Faces in outfile. Boundary coordinate increases continously along polygonal boundary. It's zero at the lexicographically smallest node on the boundary.

3853  {
3854  // Temporary storage for face elements
3855  Vector<FiniteElement*> face_el_pt;
3856 
3857  // Temporary storage for number of elements adjacent to the boundary
3858  unsigned nel = 0;
3859 
3860  // =================================================================
3861  // BEGIN: Get face elements from boundary elements
3862  // =================================================================
3863 
3864  // Temporary storage for elements adjacent to the boundary that have
3865  // an common edge (related with internal boundaries)
3866  unsigned n_repeated_ele = 0;
3867 
3868  unsigned n_regions = this->nregion();
3869 
3870 #ifdef OOMPH_HAS_MPI
3871  // map to associate the face element to the bulk element, this info.
3872  // is only necessary for the setup of boundary coordinates in a
3873  // distributed mesh when we need to extract the halo/haloed info.
3874  std::map<FiniteElement*, FiniteElement*> face_to_bulk_element_pt;
3875 #endif
3876 
3877  // Temporary storage for already done nodes
3878  Vector<std::pair<Node*, Node*>> done_nodes_pt;
3879 
3880  // If there is more than one region then only use boundary
3881  // coordinates from the bulk side (region 0)
3882  if (n_regions > 1)
3883  {
3884  for (unsigned rr = 0; rr < n_regions; rr++)
3885  {
3886  unsigned region_id = static_cast<unsigned>(this->Region_attribute[rr]);
3887 
3888 #ifdef PARANOID
3889  double diff =
3890  fabs(Region_attribute[rr] - static_cast<double>(static_cast<unsigned>(
3891  this->Region_attribute[rr])));
3892  if (diff > 0.0)
3893  {
3894  std::ostringstream error_message;
3895  error_message << "Region attributes should be unsigneds because we \n"
3896  << "only use them to set region ids\n";
3897  throw OomphLibError(error_message.str(),
3898  OOMPH_CURRENT_FUNCTION,
3899  OOMPH_EXCEPTION_LOCATION);
3900  }
3901 #endif
3902 
3903  // Loop over all elements on boundaries in region rr
3904  unsigned nel_in_region =
3905  this->nboundary_element_in_region(b, region_id);
3906  unsigned nel_repeated_in_region = 0;
3907 
3908 #ifdef PARANOID
3909  if (!Suppress_warning_about_regions_and_boundaries)
3910  {
3911  if (nel_in_region == 0)
3912  {
3913  std::ostringstream warning_message;
3914  std::string output_string =
3915  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
3916  warning_message
3917  << "There are no elements associated with boundary (" << b
3918  << ")\n"
3919  << "in region (" << region_id << "). This could happen because:\n"
3920  << "1) You did not specify boundaries with this boundary id.\n"
3921  << "---- Review carefully the indexing of your boundaries.\n"
3922  << "2) The boundary (" << b << ") is not associated with region ("
3923  << region_id << ").\n"
3924  << "---- The boundary does not touch the region.\n"
3925  << "You can suppress this warning by setting the static public "
3926  "bool\n\n"
3927  << " "
3928  "UnstructuredTwoDMeshGeometryBase::Suppress_warning_about_"
3929  "regions_and_boundaries\n\n"
3930  << "to true.\n";
3931  OomphLibWarning(
3932  warning_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
3933  }
3934  }
3935 #endif
3936 
3937  // Only bother to do anything else, if there are elements
3938  // associated with the boundary and the current region
3939  if (nel_in_region > 0)
3940  {
3941  // Flag that activates when a repeated face element is found,
3942  // possibly because we are dealing with an internal boundary
3943  bool repeated = false;
3944 
3945  // Loop over the bulk elements adjacent to boundary b
3946  for (unsigned e = 0; e < nel_in_region; e++)
3947  {
3948  // Get pointer to the bulk element that is adjacent to boundary b
3949  FiniteElement* bulk_elem_pt =
3950  this->boundary_element_in_region_pt(b, region_id, e);
3951 
3952 #ifdef OOMPH_HAS_MPI
3953  // In a distributed mesh only work with nonhalo elements
3954  if (this->is_mesh_distributed() && bulk_elem_pt->is_halo())
3955  {
3956  // Increase the number of repeated elements
3957  n_repeated_ele++;
3958  // Skip this element and go for the next one
3959  continue;
3960  }
3961 #endif
3962 
3963  // Find the index of the face of element e along boundary b
3964  int face_index =
3965  this->face_index_at_boundary_in_region(b, region_id, e);
3966 
3967  // Before adding the new element we need to be sure that
3968  // the edge that this element represent has not been
3969  // already added
3970  FiniteElement* tmp_ele_pt =
3971  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
3972 
3973  const unsigned n_nodes = tmp_ele_pt->nnode();
3974 
3975  std::pair<Node*, Node*> tmp_pair = std::make_pair(
3976  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
3977 
3978  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
3979  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
3980 
3981  // Search for repeated nodes
3982  const unsigned n_done_nodes = done_nodes_pt.size();
3983  for (unsigned l = 0; l < n_done_nodes; l++)
3984  {
3985  if (tmp_pair == done_nodes_pt[l] ||
3986  tmp_pair_inverse == done_nodes_pt[l])
3987  {
3988  nel_repeated_in_region++;
3989  repeated = true;
3990  break;
3991  }
3992  }
3993 
3994  // Create new face element?
3995  if (!repeated)
3996  {
3997  // Add the pair of nodes (edge) to the node dones
3998  done_nodes_pt.push_back(tmp_pair);
3999 #ifdef OOMPH_HAS_MPI
4000  // If the mesh is distributed then create a map from the
4001  // temporary face element to the bulk element, further
4002  // info. will be extracted from the bulk element for setup
4003  // of boundary coordinates in a distributed mesh
4004  if (this->is_mesh_distributed())
4005  {
4006  face_to_bulk_element_pt[tmp_ele_pt] = bulk_elem_pt;
4007  }
4008 #endif
4009  // Add the face element to the storage
4010  face_el_pt.push_back(tmp_ele_pt);
4011  }
4012  else
4013  {
4014  // Clean up
4015  delete tmp_ele_pt;
4016  tmp_ele_pt = 0;
4017  }
4018 
4019  // Re-start
4020  repeated = false;
4021 
4022  // Output faces?
4023  if (outfile.is_open())
4024  {
4025  face_el_pt[face_el_pt.size() - 1]->output(outfile);
4026  }
4027  } // for nel
4028 
4029  nel += nel_in_region;
4030 
4031  n_repeated_ele += nel_repeated_in_region;
4032 
4033  } // if (nel_in_region > 0)
4034 
4035  } // for (rr < n_regions)
4036 
4037  } // if (n_regions > 1)
4038  // Otherwise it's just the normal boundary functions
4039  else
4040  {
4041  // Loop over all elements on boundaries
4042  nel = this->nboundary_element(b);
4043 
4044 #ifdef PARANOID
4045  if (!Suppress_warning_about_regions_and_boundaries)
4046  {
4047  if (nel == 0)
4048  {
4049  std::ostringstream warning_message;
4050  std::string output_string =
4051  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
4052  warning_message
4053  << "There are no elements associated with boundary (" << b << ").\n"
4054  << "This could happen because you did not specify boundaries with\n"
4055  << "this boundary id. Review carefully the indexing of your\n"
4056  << "boundaries.";
4057  OomphLibWarning(
4058  warning_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4059  }
4060  }
4061 #endif
4062 
4063  // Only bother to do anything else, if there are elements
4064  if (nel > 0)
4065  {
4066  // Flag that activates when a repeated face element is found,
4067  // possibly because we are dealing with an internal boundary
4068  bool repeated = false;
4069 
4070  // Loop over the bulk elements adjacent to boundary b
4071  for (unsigned e = 0; e < nel; e++)
4072  {
4073  // Get pointer to the bulk element that is adjacent to boundary b
4074  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
4075 
4076 #ifdef OOMPH_HAS_MPI
4077 
4078  // In a distributed mesh only work with nonhalo elements
4079  if (this->is_mesh_distributed() && bulk_elem_pt->is_halo())
4080  {
4081  // Increase the number of repeated elements
4082  n_repeated_ele++;
4083  // Skip this element and go for the next one
4084  continue;
4085  }
4086 
4087 #endif
4088 
4089  // Find the index of the face of element e along boundary b
4090  int face_index = this->face_index_at_boundary(b, e);
4091 
4092  // Before adding the new element we need to be sure that the
4093  // edge that this element represent has not been already added
4094  FiniteElement* tmp_ele_pt =
4095  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
4096 
4097  const unsigned n_nodes = tmp_ele_pt->nnode();
4098 
4099  std::pair<Node*, Node*> tmp_pair = std::make_pair(
4100  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
4101 
4102  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
4103  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
4104 
4105  // Search for repeated nodes
4106  const unsigned n_done_nodes = done_nodes_pt.size();
4107  for (unsigned l = 0; l < n_done_nodes; l++)
4108  {
4109  if (tmp_pair == done_nodes_pt[l] ||
4110  tmp_pair_inverse == done_nodes_pt[l])
4111  {
4112  n_repeated_ele++;
4113  repeated = true;
4114  break;
4115  }
4116  }
4117 
4118  // Create new face element
4119  if (!repeated)
4120  {
4121  // Add the pair of nodes (edge) to the node dones
4122  done_nodes_pt.push_back(tmp_pair);
4123 #ifdef OOMPH_HAS_MPI
4124  // Create a map from the temporary face element to the bulk
4125  // element, further info. will be extracted from the bulk
4126  // element for setup of boundary coordinates in a
4127  // distributed mesh
4128  if (this->is_mesh_distributed())
4129  {
4130  face_to_bulk_element_pt[tmp_ele_pt] = bulk_elem_pt;
4131  }
4132 #endif
4133  face_el_pt.push_back(tmp_ele_pt);
4134  }
4135  else
4136  {
4137  // Free the repeated bulk element!!
4138  delete tmp_ele_pt;
4139  tmp_ele_pt = 0;
4140  }
4141 
4142  // Re-start
4143  repeated = false;
4144 
4145  // Output faces?
4146  if (outfile.is_open())
4147  {
4148  face_el_pt[face_el_pt.size() - 1]->output(outfile);
4149  }
4150 
4151  } // for (e < nel)
4152 
4153  } // if (nel > 0)
4154 
4155  } // else if (n_regions > 1)
4156 
4157  // Do not consider the repeated elements
4158  nel -= n_repeated_ele;
4159 
4160 #ifdef PARANOID
4161  if (nel != face_el_pt.size())
4162  {
4163  std::ostringstream error_message;
4164  error_message
4165  << "The independent counting of face elements (" << nel << ") for "
4166  << "boundary (" << b << ") is different\n"
4167  << "from the real number of face elements in the container ("
4168  << face_el_pt.size() << ")\n";
4169  throw OomphLibError(
4170  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4171  }
4172 #endif
4173 
4174  // =================================================================
4175  // END: Get face elements from boundary elements
4176  // =================================================================
4177 
4178  // Only bother to do anything else, if there are elements
4179  if (nel > 0)
4180  {
4181  // A flag vector to mark those face elements that are considered
4182  // as halo in the current processor
4183  std::vector<bool> is_halo_face_element(nel, false);
4184 
4185  // Count the total number of non halo face elements
4186  unsigned nnon_halo_face_elements = 0;
4187 
4188 #ifdef OOMPH_HAS_MPI
4189  // Only mark the face elements as halo if the mesh is marked as
4190  // distributed
4191  if (this->is_mesh_distributed())
4192  {
4193  for (unsigned ie = 0; ie < nel; ie++)
4194  {
4195  FiniteElement* face_ele_pt = face_el_pt[ie];
4196  // Get the bulk element
4197  FiniteElement* tmp_bulk_ele_pt = face_to_bulk_element_pt[face_ele_pt];
4198  // Check if the bulk element is halo
4199  if (!tmp_bulk_ele_pt->is_halo())
4200  {
4201  // Mark the face element as nonhalo
4202  is_halo_face_element[ie] = false;
4203  // Increase the counter for nonhalo elements
4204  nnon_halo_face_elements++;
4205  }
4206  else
4207  {
4208  // Mark the face element as halo
4209  is_halo_face_element[ie] = true;
4210  }
4211  } // for (ie < nel)
4212  } // if (this->is_mesh_distributed())
4213  else
4214  {
4215 #endif // OOMPH_HAS_MPI
4216 
4217  // If the mesh is not distributed then the number of non halo
4218  // elements is the same as the number of elements
4219  nnon_halo_face_elements = nel;
4220 
4221 #ifdef OOMPH_HAS_MPI
4222  } // else if (this->is_mesh_distributed())
4223 #endif
4224 
4225 #ifdef PARANOID
4226  // Get the total number of halo face elements
4227  const unsigned nhalo_face_element = nel - nnon_halo_face_elements;
4228 
4229  if (nhalo_face_element > 0)
4230  {
4231  std::ostringstream error_message;
4232  error_message
4233  << "There should not be halo face elements since they were not\n"
4234  << "considered when computing the face elements.\n"
4235  << "The number of found halo face elements is: ("
4236  << nhalo_face_element << ").\n\n";
4237  throw OomphLibError(error_message.str(),
4238  OOMPH_CURRENT_FUNCTION,
4239  OOMPH_EXCEPTION_LOCATION);
4240  }
4241 #endif
4242 
4243  // =================================================================
4244  // BEGIN: Sort face elements
4245  // =================================================================
4246 
4247  // The vector of lists to store the "segments" that compound the
4248  // boundaries (segments may appear only in a distributed mesh
4249  // because the boundary may have been split across multiple
4250  // processors)
4251  Vector<std::list<FiniteElement*>> segment_sorted_ele_pt;
4252 
4253  // Number of already sorted face elements (only nonhalo face
4254  // elements for a distributed mesh)
4255  unsigned nsorted_face_elements = 0;
4256 
4257  // Keep track of who's done (in a distributed mesh this apply to
4258  // nonhalo only)
4259  std::map<FiniteElement*, bool> done_el;
4260 
4261  // Keep track of which element is inverted (in distributed mesh
4262  // the elements may be inverted with respect to the segment they
4263  // belong)
4264  std::map<FiniteElement*, bool> is_inverted;
4265 
4266  // Iterate until all possible segments have been created. In a
4267  // non distributed mesh there is only one segment which defines
4268  // the complete boundary
4269  while (nsorted_face_elements < nnon_halo_face_elements)
4270  {
4271  // The sorted list of face elements (in a distributed mesh a
4272  // collection of continuous face elements define a segment)
4273  std::list<FiniteElement*> sorted_el_pt;
4274 
4275  FiniteElement* ele_face_pt = 0;
4276 
4277 #ifdef PARANOID
4278  // Select an initial element for the segment
4279  bool found_initial_face_element = false;
4280 #endif
4281 
4282  // Store the index of the initial face element
4283  unsigned iface = 0;
4284 #ifdef OOMPH_HAS_MPI
4285  if (this->is_mesh_distributed())
4286  {
4287  for (iface = 0; iface < nel; iface++)
4288  {
4289  // Only work with nonhalo face elements
4290  if (!is_halo_face_element[iface])
4291  {
4292  ele_face_pt = face_el_pt[iface];
4293  // If not done then take it as initial face element
4294  if (!done_el[ele_face_pt])
4295  {
4296 #ifdef PARANOID
4297  // Set the flag to indicate the initial element was
4298  // found
4299  found_initial_face_element = true;
4300 #endif
4301  // Increase the number of sorted face elements
4302  nsorted_face_elements++;
4303  // Set the index to the next face element
4304  iface++;
4305  // Add the face element in the container
4306  sorted_el_pt.push_back(ele_face_pt);
4307  // Mark as done
4308  done_el[ele_face_pt] = true;
4309  break;
4310  } // if (!done_el[ele_face_pt])
4311  } // if (!is_halo_face_element[iface])
4312  } // for (iface < nel)
4313  } // if (this->is_mesh_distributed())
4314  else
4315  {
4316 #endif // #ifdef OOMPH_HAS_MPI
4317 
4318  // When the mesh is not distributed just take the first
4319  // element and put it in the sorted list
4320  ele_face_pt = face_el_pt[0];
4321 #ifdef PARANOID
4322  // Set the flag to indicate the initial element was found
4323  found_initial_face_element = true;
4324 #endif
4325  // Increase the number of sorted face elements
4326  nsorted_face_elements++;
4327  // Set the index to the next face element
4328  iface = 1;
4329  // Add the face element in the container
4330  sorted_el_pt.push_back(ele_face_pt);
4331  // Mark as done
4332  done_el[ele_face_pt] = true;
4333 
4334 #ifdef OOMPH_HAS_MPI
4335  } // else if (this->is_mesh_distributed())
4336 #endif
4337 
4338 #ifdef PARANOID
4339  if (!found_initial_face_element)
4340  {
4341  std::ostringstream error_message;
4342  error_message << "Could not find an initial face element for the "
4343  "current segment\n";
4344  throw OomphLibError(error_message.str(),
4345  OOMPH_CURRENT_FUNCTION,
4346  OOMPH_EXCEPTION_LOCATION);
4347  }
4348 #endif
4349 
4350  // Number of nodes of the initial face element
4351  const unsigned nnod = ele_face_pt->nnode();
4352 
4353  // Left and rightmost nodes (the left and right nodes of the
4354  // current face element)
4355  Node* left_node_pt = ele_face_pt->node_pt(0);
4356  Node* right_node_pt = ele_face_pt->node_pt(nnod - 1);
4357 
4358  // Continue iterating if a new face element has been added to
4359  // the list
4360  bool face_element_added = false;
4361 
4362  // While a new face element has been added to the set of sorted
4363  // face elements continue iterating
4364  do
4365  {
4366  // Start from the next face element since we have already
4367  // added the previous one as the initial face element (any
4368  // previous face element had to be added on previous
4369  // iterations)
4370  for (unsigned iiface = iface; iiface < nel; iiface++)
4371  {
4372  // Re-start flag
4373  face_element_added = false;
4374 
4375  // Get the candidate element
4376  ele_face_pt = face_el_pt[iiface];
4377 
4378  // Check that the candidate element has not been done and
4379  // is not a halo element
4380  if (!(done_el[ele_face_pt] || is_halo_face_element[iiface]))
4381  {
4382  // Get the left and right nodes of the current element
4383  Node* local_left_node_pt = ele_face_pt->node_pt(0);
4384  Node* local_right_node_pt = ele_face_pt->node_pt(nnod - 1);
4385 
4386  // New element fits at the left of segment and is not inverted
4387  if (left_node_pt == local_right_node_pt)
4388  {
4389  left_node_pt = local_left_node_pt;
4390  sorted_el_pt.push_front(ele_face_pt);
4391  is_inverted[ele_face_pt] = false;
4392  face_element_added = true;
4393  }
4394  // New element fits at the left of segment and is inverted
4395  else if (left_node_pt == local_left_node_pt)
4396  {
4397  left_node_pt = local_right_node_pt;
4398  sorted_el_pt.push_front(ele_face_pt);
4399  is_inverted[ele_face_pt] = true;
4400  face_element_added = true;
4401  }
4402  // New element fits on the right of segment and is not inverted
4403  else if (right_node_pt == local_left_node_pt)
4404  {
4405  right_node_pt = local_right_node_pt;
4406  sorted_el_pt.push_back(ele_face_pt);
4407  is_inverted[ele_face_pt] = false;
4408  face_element_added = true;
4409  }
4410  // New element fits on the right of segment and is inverted
4411  else if (right_node_pt == local_right_node_pt)
4412  {
4413  right_node_pt = local_left_node_pt;
4414  sorted_el_pt.push_back(ele_face_pt);
4415  is_inverted[ele_face_pt] = true;
4416  face_element_added = true;
4417  }
4418 
4419  if (face_element_added)
4420  {
4421  done_el[ele_face_pt] = true;
4422  nsorted_face_elements++;
4423  break;
4424  }
4425 
4426  } // if (!(done_el[ele_face_pt] || is_halo_face_element[iiface]))
4427  } // for (iiface<nnon_halo_face_element)
4428  } while (face_element_added &&
4429  (nsorted_face_elements < nnon_halo_face_elements));
4430 
4431  // Store the created segment in the vector of segments
4432  segment_sorted_ele_pt.push_back(sorted_el_pt);
4433 
4434  } // while(nsorted_face_elements < nnon_halo_face_elements);
4435 
4436 #ifdef OOMPH_HAS_MPI
4437  if (!this->is_mesh_distributed())
4438  {
4439 #endif
4440  // Are we done?
4441  if (nsorted_face_elements != nel || segment_sorted_ele_pt.size() != 1)
4442  {
4443  std::ostringstream error_message;
4444  error_message << "Was only able to setup boundary coordinate on "
4445  << "boundary " << b << "\nfor " << nsorted_face_elements
4446  << " of " << nel
4447  << " face elements. This usually means\n"
4448  << "that the boundary is not simply connected.\n\n"
4449  << "Re-run the setup_boundary_coordintes() function\n"
4450  << "with an output file specified "
4451  << "as the second argument.\n"
4452  << "This file will contain FaceElements that\n"
4453  << "oomph-lib believes to be located on the boundary.\n"
4454  << std::endl;
4455  throw OomphLibError(error_message.str(),
4456  OOMPH_CURRENT_FUNCTION,
4457  OOMPH_EXCEPTION_LOCATION);
4458  }
4459 #ifdef OOMPH_HAS_MPI
4460  } // if (!this->is_mesh_distributed())
4461 #endif
4462 
4463  // =================================================================
4464  // END: Sort face elements
4465  // =================================================================
4466 
4467  // ----------------------------------------------------------------
4468 
4469  // =================================================================
4470  // BEGIN: Assign global/local (non distributed mesh/distributed
4471  // mesh) boundary coordinates to nodes
4472  // =================================================================
4473 
4474  // Compute the (local) boundary coordinates of the nodes in the
4475  // segments. In a distributed mesh this info. will be used by a
4476  // root processor to compute the (global) boundary coordinates.
4477 
4478  // Vector of sets that stores the nodes of each segment based on
4479  // a lexicographically order starting from the bottom left node
4480  // of each segment
4481  Vector<std::set<Node*>> segment_all_nodes_pt;
4482 
4483  // The number of segments in this processor
4484  const unsigned nsegments = segment_sorted_ele_pt.size();
4485 
4486 #ifdef PARANOID
4487  if (nnon_halo_face_elements > 0 && nsegments == 0)
4488  {
4489  std::ostringstream error_message;
4490  error_message
4491  << "The number of segments is zero, but the number of nonhalo\n"
4492  << "elements is: (" << nnon_halo_face_elements << ")\n";
4493  throw OomphLibError(error_message.str(),
4494  OOMPH_CURRENT_FUNCTION,
4495  OOMPH_EXCEPTION_LOCATION);
4496  } // if (nnon_halo_face_elements > 0 && nsegments == 0)
4497 
4498 #endif
4499 
4500  // Store the arclength of each segment in the current processor
4501  Vector<double> segment_arclength(nsegments);
4502 
4503 #ifdef OOMPH_HAS_MPI
4504 
4505  // Clear the boundary segment nodes storage
4506  this->flush_boundary_segment_node(b);
4507 
4508  // Set the number of segments for the current boundary
4509  this->set_nboundary_segment_node(b, nsegments);
4510 
4511 #endif // #ifdef OOMPH_HAS_MPI
4512 
4513  // Go through all the segments and compute their (local) boundary
4514  // coordinates
4515  for (unsigned is = 0; is < nsegments; is++)
4516  {
4517 #ifdef PARANOID
4518 
4519  if (segment_sorted_ele_pt[is].size() == 0)
4520  {
4521  std::ostringstream error_message;
4522  std::string output_string =
4523  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
4524  error_message << "The (" << is << ")-th segment has no elements\n";
4525  throw OomphLibError(
4526  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4527  } // if (segment_sorted_ele_pt[is].size() == 0)
4528 
4529 #endif
4530 
4531  // Get access to the first element on the segment
4532  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
4533 
4534  // Number of nodes
4535  unsigned nnod = first_ele_pt->nnode();
4536 
4537  // Get the first node of the current segment
4538  Node* first_node_pt = first_ele_pt->node_pt(0);
4539  if (is_inverted[first_ele_pt])
4540  {
4541  first_node_pt = first_ele_pt->node_pt(nnod - 1);
4542  }
4543 
4544  // Coordinates of left node
4545  double x_left = first_node_pt->x(0);
4546  double y_left = first_node_pt->x(1);
4547 
4548  // Initialise boundary coordinate (local boundary coordinate
4549  // for boundaries with more than one segment)
4550  Vector<double> zeta(1, 0.0);
4551 
4552  // Set boundary coordinate
4553  first_node_pt->set_coordinates_on_boundary(b, zeta);
4554 
4555  // Lexicographically bottom left node
4556  std::set<Node*> local_nodes_pt;
4557  local_nodes_pt.insert(first_node_pt);
4558 
4559 #ifdef OOMPH_HAS_MPI
4560 
4561  // Insert the node in the look-up scheme for nodes in segments
4562  this->add_boundary_segment_node(b, is, first_node_pt);
4563 
4564 #endif // #ifdef OOMPH_HAS_MPI
4565 
4566  // Now loop over nodes in order
4567  for (std::list<FiniteElement*>::iterator it =
4568  segment_sorted_ele_pt[is].begin();
4569  it != segment_sorted_ele_pt[is].end();
4570  it++)
4571  {
4572  // Get element
4573  FiniteElement* el_pt = *it;
4574 
4575  // Start node and increment
4576  unsigned k_nod = 1;
4577  int nod_diff = 1;
4578  if (is_inverted[el_pt])
4579  {
4580  k_nod = nnod - 2;
4581  nod_diff = -1;
4582  }
4583 
4584  // Loop over nodes
4585  for (unsigned j = 1; j < nnod; j++)
4586  {
4587  Node* nod_pt = el_pt->node_pt(k_nod);
4588  k_nod += nod_diff;
4589 
4590  // Coordinates of right node
4591  double x_right = nod_pt->x(0);
4592  double y_right = nod_pt->x(1);
4593 
4594  // Increment boundary coordinate
4595  zeta[0] += sqrt((x_right - x_left) * (x_right - x_left) +
4596  (y_right - y_left) * (y_right - y_left));
4597 
4598  // Set boundary coordinate
4599  nod_pt->set_coordinates_on_boundary(b, zeta);
4600 
4601  // Increment reference coordinate
4602  x_left = x_right;
4603  y_left = y_right;
4604 
4605  // Get lexicographically bottom left node but only
4606  // use vertex nodes as candidates
4607  local_nodes_pt.insert(nod_pt);
4608 
4609 #ifdef OOMPH_HAS_MPI
4610 
4611  // Insert the node in the look-up scheme for nodes in segments
4612  this->add_boundary_segment_node(b, is, nod_pt);
4613 
4614 #endif // #ifdef OOMPH_HAS_MPI
4615 
4616  } // for (j < nnod)
4617  } // iterator over the elements in the segment
4618 
4619  // Assign the arclength of the current segment
4620  segment_arclength[is] = zeta[0];
4621 
4622  // Add the nodes for the corresponding segment in the container
4623  segment_all_nodes_pt.push_back(local_nodes_pt);
4624 
4625  } // for (is < nsegments)
4626 
4627  // =================================================================
4628  // END: Assign global/local (non distributed mesh/distributed
4629  // mesh) boundary coordinates to nodes
4630  // =================================================================
4631 
4632  // Store the initial arclength for each segment of boundary in
4633  // the current processor, initialise to zero in case we have a
4634  // non distributed mesh
4635  Vector<double> initial_segment_arclength(nsegments, 0.0);
4636 
4637  // Store the final arclength for each segment of boundary in the
4638  // current processor, initalise to zero in case we have a non
4639  // distributed mesh
4640  Vector<double> final_segment_arclength(nsegments, 0.0);
4641 
4642  // Store the initial zeta for each segment of boundary in the
4643  // current processor, initalise to zero in case we have a non
4644  // distributed mesh
4645  Vector<double> initial_segment_zeta(nsegments, 0.0);
4646 
4647  // Store the final zeta for each segment of boundary in the
4648  // current processor, initalise to zero in case we have a non
4649  // distributed mesh
4650  Vector<double> final_segment_zeta(nsegments, 0.0);
4651 
4652  // --------------------------------------------------------------
4653  // DISTRIBUTED MESH: BEGIN - Check orientation of boundaries and
4654  // assign boundary coordinates accordingly
4655  // --------------------------------------------------------------
4656 
4657 #ifdef OOMPH_HAS_MPI
4658 
4659  // Check that the mesh is distributed and that the initial zeta
4660  // values for the boundary segments have been already
4661  // assigned. When the method is called by the first time (when
4662  // the mesh is just created) the initial zeta values for the
4663  // boundary segments are not available
4664  if (this->is_mesh_distributed() &&
4665  Assigned_segments_initial_zeta_values[b])
4666  {
4667  // Get the initial and final coordinates of each segment
4668 
4669  // For each segment in the processor check whether it was
4670  // inverted in the root processor, if that is the case then it
4671  // is necessary to re-sort the face elements and invert them
4672  for (unsigned is = 0; is < nsegments; is++)
4673  {
4674  // Check if we need/or not to invert the current zeta values
4675  // on the segment (only apply for boundaries with GeomObject
4676  // associated)
4677 
4678  // Does the boundary HAS a GeomObject associated?
4679  if (this->boundary_geom_object_pt(b) != 0)
4680  {
4681  // Get the first and last node of the current segment and
4682  // their zeta values
4683 
4684  // Get access to the first element on the segment
4685  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
4686 
4687  // Number of nodes
4688  const unsigned nnod = first_ele_pt->nnode();
4689 
4690  // Get the first node of the current segment
4691  Node* first_node_pt = first_ele_pt->node_pt(0);
4692  if (is_inverted[first_ele_pt])
4693  {
4694  first_node_pt = first_ele_pt->node_pt(nnod - 1);
4695  }
4696 
4697  // Get access to the last element on the segment
4698  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
4699 
4700  // Get the last node of the current segment
4701  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
4702  if (is_inverted[last_ele_pt])
4703  {
4704  last_node_pt = last_ele_pt->node_pt(0);
4705  }
4706 
4707  // Get the zeta coordinates for the first and last node
4708  Vector<double> current_segment_initial_zeta(1);
4709  Vector<double> current_segment_final_zeta(1);
4710 
4711  first_node_pt->get_coordinates_on_boundary(
4712  b, current_segment_initial_zeta);
4713  last_node_pt->get_coordinates_on_boundary(
4714  b, current_segment_final_zeta);
4715 
4716 #ifdef PARANOID
4717  // Check whether the zeta values in the segment are in
4718  // increasing or decreasing order
4719  if (current_segment_initial_zeta[0] < current_segment_final_zeta[0])
4720  {
4721  // do nothing
4722  }
4723  else if (current_segment_initial_zeta[0] >
4724  current_segment_final_zeta[0])
4725  {
4726  std::stringstream error_message;
4727  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
4728  "setup_boundary_coordinates()";
4729  error_message
4730  << "The zeta values are in decreasing order, this is weird\n"
4731  << "since they have just been set-up in increasing order few\n"
4732  << "lines above\n"
4733  << "Boundary: (" << b << ")\n"
4734  << "Segment: (" << is << ")\n"
4735  << "Initial zeta value: (" << current_segment_initial_zeta[0]
4736  << ")\n"
4737  << "Initial coordinate: (" << first_node_pt->x(0) << ", "
4738  << first_node_pt->x(1) << ")\n"
4739  << "Final zeta value: (" << current_segment_final_zeta[0]
4740  << ")\n"
4741  << "Initial coordinate: (" << last_node_pt->x(0) << ", "
4742  << last_node_pt->x(1) << ")\n";
4743  throw OomphLibError(
4744  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4745  }
4746  else
4747  {
4748  std::stringstream error_message;
4749  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
4750  "setup_boundary_coordinates()";
4751  error_message
4752  << "It was not possible to determine whether the zeta values on"
4753  << "the current segment\nof the boundary are in"
4754  << "increasing or decreasing order\n\n"
4755  << "Boundary: (" << b << ")\n"
4756  << "Segment: (" << is << ")\n"
4757  << "Arclength: (" << segment_arclength[is] << "\n"
4758  << "Initial zeta value: (" << current_segment_initial_zeta[0]
4759  << ")\n"
4760  << "Initial coordinate: (" << first_node_pt->x(0) << ", "
4761  << first_node_pt->x(1) << ")\n"
4762  << "Final zeta value: (" << current_segment_final_zeta[0]
4763  << ")\n"
4764  << "Initial coordinate: (" << last_node_pt->x(0) << ", "
4765  << last_node_pt->x(1) << ")\n";
4766  throw OomphLibError(
4767  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4768  }
4769 #endif // #ifdef PARANOID
4770 
4771  // Now get the original zeta values and check if they are in
4772  // increasing or decreasing order
4773  const double original_segment_initial_zeta =
4774  boundary_segment_initial_zeta(b)[is];
4775  const double original_segment_final_zeta =
4776  boundary_segment_final_zeta(b)[is];
4777 
4778  // Now check if the zeta values go in increase or decrease
4779  // order
4780  if (original_segment_final_zeta > original_segment_initial_zeta)
4781  {
4782  // The original values go in increasing order, only need
4783  // to change the values if the original segment is marked
4784  // as inverted
4785  if (boundary_segment_inverted(b)[is])
4786  {
4787  // The original segment is inverted, then we need to
4788  // reverse the boundary coordinates. Go through all the
4789  // nodes and reverse its value
4790  std::set<Node*> all_nodes_pt = segment_all_nodes_pt[is];
4791  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
4792  it != all_nodes_pt.end();
4793  it++)
4794  {
4795  Vector<double> zeta(1);
4796  // Get the node
4797  Node* nod_pt = (*it);
4798  // Get the boundary coordinate associated to the node
4799  nod_pt->get_coordinates_on_boundary(b, zeta);
4800  // Compute its new value
4801  zeta[0] = segment_arclength[is] - zeta[0];
4802  // Set the new boundary coordinate value
4803  nod_pt->set_coordinates_on_boundary(b, zeta);
4804  } // Loop over the nodes in the segment
4805  } // if (boundary_segment_inverted(b)[is])
4806  else
4807  {
4808  // The boundary is not inverted, do nothing!!!
4809  }
4810 
4811  } // original zeta values in increasing order
4812  else if (original_segment_final_zeta <
4813  original_segment_initial_zeta)
4814  {
4815  // The original values go in decreasing order, only need
4816  // to change the values if the original segment is NOT
4817  // marked as inverted
4818  if (boundary_segment_inverted(b)[is])
4819  {
4820  // The boundary is inverted, do nothing!!!
4821  }
4822  else
4823  {
4824  // The original segment is NOT inverted, then we need
4825  // to reverse the boundary coordinates. Go through all
4826  // the nodes and reverse its value
4827  std::set<Node*> all_nodes_pt = segment_all_nodes_pt[is];
4828  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
4829  it != all_nodes_pt.end();
4830  it++)
4831  {
4832  Vector<double> zeta(1);
4833  // Get the node
4834  Node* nod_pt = (*it);
4835  // Get the boundary coordinate associated to the node
4836  nod_pt->get_coordinates_on_boundary(b, zeta);
4837  // Compute its new value
4838  zeta[0] = segment_arclength[is] - zeta[0];
4839  // Set the new boundary coordinate value
4840  nod_pt->set_coordinates_on_boundary(b, zeta);
4841  } // Loop over the nodes in the segment
4842  } // else if (boundary_segment_inverted(b)[is])
4843  } // original zeta values in decreasing order
4844 #ifdef PARANOID
4845  else
4846  {
4847  std::stringstream error_message;
4848  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
4849  "setup_boundary_coordinates()";
4850  error_message
4851  << "It was not possible to identify if the zeta values on the\n"
4852  << "current segment in the boundary should go in increasing\n"
4853  << "or decreasing order.\n\n"
4854  << "Boundary: (" << b << ")\n"
4855  << "Segment: (" << is << ")\n"
4856  << "Initial zeta value: (" << original_segment_initial_zeta
4857  << ")\n"
4858  << "Final zeta value: (" << original_segment_final_zeta
4859  << ")\n";
4860  throw OomphLibError(
4861  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4862  }
4863 #endif
4864 
4865  } // if (this->boundary_geom_object_pt(b)!=0)
4866  else
4867  {
4868  // No GeomObject associated, do nothing!!!
4869 
4870  } // else if (this->boundary_geom_object_pt(b)!=0)
4871  } // for (is < nsegments)
4872  } // if (this->is_mesh_distributed() &&
4873  // Assigned_segments_initial_zeta_values[b])
4874 
4875 #endif // #ifdef OOMPH_HAS_MPI
4876 
4877  // Get the number of sets for nodes
4878 #ifdef PARANOID
4879  if (segment_all_nodes_pt.size() != nsegments)
4880  {
4881  std::ostringstream error_message;
4882  std::string output_string =
4883  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
4884  error_message << "The number of segments (" << nsegments
4885  << ") and the number of "
4886  << "sets of nodes (" << segment_all_nodes_pt.size()
4887  << ") representing\n"
4888  << "the\nsegments is different!!!\n\n";
4889  throw OomphLibError(
4890  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4891  }
4892 #endif
4893 
4894  // --------------------------------------------------------------
4895  // DISTRIBUTED MESH: END - Check orientation of boundaries and
4896  // assign boundary coordinates accordingly
4897  // --------------------------------------------------------------
4898 
4899  // =================================================================
4900  // BEGIN: Get the lenght of the boundaries or segments of the
4901  // boundary
4902  // =================================================================
4903 
4904  // The nodes have been assigned arc-length coordinates from one
4905  // end or the other of the segments. If the mesh is distributed
4906  // the values are set so that they agree with the increasing or
4907  // decreasing order of the zeta values for the segments
4908 
4909  // Storage for the coordinates of the first and last nodes
4910  Vector<double> first_coordinate(2);
4911  Vector<double> last_coordinate(2);
4912  // Storage for the zeta coordinate of the first and last nodes
4913  Vector<double> first_node_zeta_coordinate(1, 0.0);
4914  Vector<double> last_node_zeta_coordinate(1, 0.0);
4915 
4916  // Store the accumulated arclength, used at the end to check the
4917  // max arclength of the boundary
4918  double boundary_arclength = 0.0;
4919 
4920  // If the mesh is marked as distributed and the initial zeta
4921  // values for the segments have been computed then get the
4922  // info. regarding the initial and final nodes coordinates, same
4923  // as the zeta boundary values for those nodes
4924 
4925 #ifdef OOMPH_HAS_MPI
4926  if (this->is_mesh_distributed() &&
4927  Assigned_segments_initial_zeta_values[b])
4928  {
4929  // --------------------------------------------------------------
4930  // DISTRIBUTED MESH: BEGIN
4931  // --------------------------------------------------------------
4932 
4933  // Get the initial and final coordinates for the complete boundary
4934  first_coordinate = boundary_initial_coordinate(b);
4935  last_coordinate = boundary_final_coordinate(b);
4936  // Get the initial and final zeta values for the boundary
4937  // (arclength)
4938  first_node_zeta_coordinate = boundary_initial_zeta_coordinate(b);
4939  last_node_zeta_coordinate = boundary_final_zeta_coordinate(b);
4940 
4941  // The total arclength is the maximum between the initial and
4942  // final zeta coordinate
4943  boundary_arclength =
4944  std::max(first_node_zeta_coordinate[0], last_node_zeta_coordinate[0]);
4945 
4946 #ifdef PARANOID
4947  if (boundary_arclength == 0)
4948  {
4949  std::ostringstream error_message;
4950  std::string output_string =
4951  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
4952  error_message << "The boundary arclength is zero for boundary (" << b
4953  << ")\n";
4954  throw OomphLibError(
4955  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4956  }
4957 #endif
4958 
4959  // Check if there is a GeomObject associated to the boundary,
4960  // and assign the corresponding zeta (arclength) values
4961  if (this->boundary_geom_object_pt(b) != 0)
4962  {
4963  initial_segment_zeta = boundary_segment_initial_zeta(b);
4964  final_segment_zeta = boundary_segment_final_zeta(b);
4965  }
4966  else
4967  {
4968  initial_segment_arclength = boundary_segment_initial_arclength(b);
4969  final_segment_arclength = boundary_segment_final_arclength(b);
4970  }
4971 
4972  // --------------------------------------------------------------
4973  // DISTRIBUTED MESH: END
4974  // --------------------------------------------------------------
4975 
4976  } // if (this->is_mesh_distributed() &&
4977  // Assigned_segments_initial_zeta_values[b])
4978  else
4979  {
4980 #endif // #ifdef OOMPH_HAS_MPI
4981 
4982  // If the mesh is NOT distributed or the initial and final zeta
4983  // values of the segments have not been assigned then perform
4984  // as in the serial case
4985  if (this->boundary_geom_object_pt(b) != 0)
4986  {
4987  initial_segment_zeta[0] = this->boundary_coordinate_limits(b)[0];
4988  final_segment_zeta[0] = this->boundary_coordinate_limits(b)[1];
4989  }
4990  else
4991  {
4992  // When the mesh is or not distributed but the initial
4993  // segment's zeta values HAVE NOT been established then
4994  // initalize the initial segment to zero
4995  initial_segment_arclength[0] = 0.0;
4996  }
4997 #ifdef OOMPH_HAS_MPI
4998  } // The mesh is NOT distributed or the zeta values for the
4999  // segments have not been established
5000 
5001 #endif // #ifdef OOMPH_HAS_MPI
5002 
5003  // =================================================================
5004  // END: Get the lenght of the boundaries or segments of the
5005  // boundary
5006  // =================================================================
5007 
5008  // =================================================================
5009  // BEGIN: Scale the boundary coordinates based on whether the
5010  // boundary has an associated GeomObj or not
5011  // =================================================================
5012 
5013  // Go through all the segments and assign the scaled boundary
5014  // coordinates
5015  for (unsigned is = 0; is < nsegments; is++)
5016  {
5017 #ifdef PARANOID
5018  if (segment_sorted_ele_pt[is].size() == 0)
5019  {
5020  std::ostringstream error_message;
5021  std::string output_string =
5022  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
5023  error_message << "The (" << is << ")-th segment has no elements\n";
5024  throw OomphLibError(
5025  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
5026  } // if (segment_sorted_ele_pt[is].size() == 0)x
5027 #endif
5028 
5029  // Get the first and last nodes coordinates for the current
5030  // segment
5031 
5032 #ifdef OOMPH_HAS_MPI
5033  // Check if the initial zeta values of the segments have been
5034  // established, if they have not then get the first and last
5035  // coordinates from the current segments, same as the zeta
5036  // values
5037  if (!Assigned_segments_initial_zeta_values[b])
5038  {
5039 #endif // #ifdef OOMPH_HAS_MPI
5040 
5041  // Get access to the first element on the segment
5042  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
5043 
5044  // Number of nodes
5045  const unsigned nnod = first_ele_pt->nnode();
5046 
5047  // Get the first node of the current segment
5048  Node* first_node_pt = first_ele_pt->node_pt(0);
5049  if (is_inverted[first_ele_pt])
5050  {
5051  first_node_pt = first_ele_pt->node_pt(nnod - 1);
5052  }
5053 
5054  // Get access to the last element on the segment
5055  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
5056 
5057  // Get the last node of the current segment
5058  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
5059  if (is_inverted[last_ele_pt])
5060  {
5061  last_node_pt = last_ele_pt->node_pt(0);
5062  }
5063 
5064  // Get the coordinates for the first and last node
5065  for (unsigned i = 0; i < 2; i++)
5066  {
5067  first_coordinate[i] = first_node_pt->x(i);
5068  last_coordinate[i] = last_node_pt->x(i);
5069  }
5070 
5071  // Get the zeta coordinates for the first and last node
5072  first_node_pt->get_coordinates_on_boundary(
5073  b, first_node_zeta_coordinate);
5074  last_node_pt->get_coordinates_on_boundary(b,
5075  last_node_zeta_coordinate);
5076 
5077 #ifdef OOMPH_HAS_MPI
5078  } // if (!Assigned_segments_initial_zeta_values[b])
5079 #endif // #ifdef OOMPH_HAS_MPI
5080 
5081  // Get the nodes of the current segment
5082  std::set<Node*> all_nodes_pt = segment_all_nodes_pt[is];
5083 
5084  // If the boundary has a geometric object representation then
5085  // scale the coordinates to match those of the geometric object
5086  GeomObject* const geom_object_pt = this->boundary_geom_object_pt(b);
5087  if (geom_object_pt != 0)
5088  {
5089  Vector<double> bound_coord_limits =
5090  this->boundary_coordinate_limits(b);
5091  // Get the position of the ends of the geometric object
5092  Vector<double> zeta(1);
5093  Vector<double> first_geom_object_location(2);
5094  Vector<double> last_geom_object_location(2);
5095 
5096  // Get the zeta value for the initial coordinates of the
5097  // GeomObject
5098  zeta[0] = bound_coord_limits[0];
5099  // zeta[0] = initial_segment_zeta[is];
5100 
5101  // Get the coordinates from the initial zeta value from the
5102  // GeomObject
5103  geom_object_pt->position(zeta, first_geom_object_location);
5104 
5105  // Get the zeta value for the final coordinates of the
5106  // GeomObject
5107  zeta[0] = bound_coord_limits[1];
5108  // zeta[0] = final_segment_zeta[is];
5109 
5110  // Get the coordinates from the final zeta value from the
5111  // GeomObject
5112  geom_object_pt->position(zeta, last_geom_object_location);
5113 
5114  // Calculate the errors in position between the first and
5115  // last nodes and the endpoints of the geometric object
5116  double error = 0.0;
5117  double tmp_error = 0.0;
5118  for (unsigned i = 0; i < 2; i++)
5119  {
5120  const double dist =
5121  first_geom_object_location[i] - first_coordinate[i];
5122  tmp_error += dist * dist;
5123  }
5124  error += sqrt(tmp_error);
5125  tmp_error = 0.0;
5126  for (unsigned i = 0; i < 2; i++)
5127  {
5128  const double dist =
5129  last_geom_object_location[i] - last_coordinate[i];
5130  tmp_error += dist * dist;
5131  }
5132  error += sqrt(tmp_error);
5133 
5134  // Calculate the errors in position between the first and
5135  // last nodes and the endpoints of the geometric object if
5136  // reversed
5137  double rev_error = 0.0;
5138  tmp_error = 0.0;
5139  for (unsigned i = 0; i < 2; i++)
5140  {
5141  const double dist =
5142  first_geom_object_location[i] - last_coordinate[i];
5143  tmp_error += dist * dist;
5144  }
5145  rev_error += sqrt(tmp_error);
5146  tmp_error = 0.0;
5147  for (unsigned i = 0; i < 2; i++)
5148  {
5149  const double dist =
5150  last_geom_object_location[i] - first_coordinate[i];
5151  tmp_error += dist * dist;
5152  }
5153  rev_error += sqrt(tmp_error);
5154 
5155  // Number of polyline vertices along this boundary
5156  const unsigned n_vertex = Polygonal_vertex_arclength_info[b].size();
5157 
5158  // Get polygonal vertex data
5159  Vector<std::pair<double, double>> polygonal_vertex_arclength =
5160  Polygonal_vertex_arclength_info[b];
5161 
5162  // If the (normal) error is small than reversed then we have
5163  // the coordinate direction correct. If not then we must
5164  // reverse it bool reversed = false;
5165  if (error < rev_error)
5166  {
5167  // Coordinates are aligned (btw: don't delete this block --
5168  // there's a final else below to catch errors!)
5169 
5170  // reversed = false;
5171  }
5172  else if (error > rev_error)
5173  {
5174  // reversed = true;
5175 
5176  // Reverse the limits of the boundary coordinates along the
5177  // geometric object
5178  double temp = bound_coord_limits[0];
5179  bound_coord_limits[0] = bound_coord_limits[1];
5180  bound_coord_limits[1] = temp;
5181 
5182 #ifdef OOMPH_HAS_MPI
5183  // If we are working with a NON distributed mesh then
5184  // re-assign the initial and final zeta values
5185  if (!this->is_mesh_distributed())
5186  {
5187 #endif
5188  temp = initial_segment_zeta[is];
5189  initial_segment_zeta[is] = final_segment_zeta[is];
5190  final_segment_zeta[is] = temp;
5191 #ifdef OOMPH_HAS_MPI
5192  }
5193 #endif
5194  // Reverse the vertices information
5195  for (unsigned v = 0; v < n_vertex; v++)
5196  {
5197  polygonal_vertex_arclength[v].first =
5198  Polygonal_vertex_arclength_info[b][v].first;
5199 
5200  polygonal_vertex_arclength[v].second =
5201  Polygonal_vertex_arclength_info[b][n_vertex - v - 1].second;
5202  } // for (v < n_vertex)
5203  }
5204  else
5205  {
5206  std::ostringstream error_stream;
5207  std::string output_string =
5208  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
5209  error_stream
5210  << "Something very strange has happened.\n"
5211  << "The error between the endpoints of the geometric object\n"
5212  << "and the first and last nodes on the boundary is the same\n"
5213  << "irrespective of the direction of the coordinate.\n"
5214  << "This probably means that things are way off.\n"
5215  << "The errors are " << error << " and " << rev_error << "\n";
5216  std::cout << error_stream.str();
5217  throw OomphLibError(
5218  error_stream.str(), output_string, OOMPH_EXCEPTION_LOCATION);
5219  }
5220 
5221  // Get the total arclength of the edge
5222  // last_node_pt->get_coordinates_on_boundary(b, zeta);
5223  // double zeta_old_range=zeta[0];
5224 
5225  // Store the arclength of the segment (not yet mapped to the
5226  // boundary coordinates defined by the GeomObject)
5227  const double zeta_old_range = segment_arclength[is];
5228 
5229  // double zeta_new_range=bound_coord_limits[1]-bound_coord_limits[0];
5230 
5231  // The range to map the zeta values
5232  double zeta_new_range =
5233  final_segment_zeta[is] - initial_segment_zeta[is];
5234  // The initial zeta value for the current segment
5235  double initial_local_segment_zeta = initial_segment_zeta[is];
5236 
5237 #ifdef OOMPH_HAS_MPI
5238 
5239  // If we are working with a distributed mes and the initial
5240  // and final zeta values for the segments have been
5241  // established then reset the range to map
5242  if (this->is_mesh_distributed() &&
5243  Assigned_segments_initial_zeta_values[b])
5244  {
5245  // Re-set the range to map the zeta values
5246  zeta_new_range =
5247  std::fabs(final_segment_zeta[is] - initial_segment_zeta[is]);
5248 
5249  // Re-set the initial zeta value of the segment
5250  initial_local_segment_zeta =
5251  std::min(initial_segment_zeta[is], final_segment_zeta[is]);
5252  }
5253 
5254 #endif
5255 
5256  // Re-assign boundary coordinate for the case where boundary
5257  // is represented by polygon
5258  unsigned use_old = false;
5259  if (n_vertex == 0) use_old = true;
5260 
5261  // Now scale the coordinates accordingly
5262  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
5263  it != all_nodes_pt.end();
5264  it++)
5265  {
5266  // Get the node
5267  Node* nod_pt = (*it);
5268 
5269  // Get coordinate based on arclength along polygonal repesentation
5270  nod_pt->get_coordinates_on_boundary(b, zeta);
5271 
5272  if (use_old)
5273  {
5274  // Boundary is actually a polygon -- simply rescale
5275  zeta[0] = initial_local_segment_zeta +
5276  (zeta_new_range / zeta_old_range) * (zeta[0]);
5277  }
5278  else
5279  {
5280  // Scale such that vertex nodes stay where they were
5281  bool found = false;
5282 
5283  // Loop over vertex nodes
5284  for (unsigned v = 1; v < n_vertex; v++)
5285  {
5286  if ((zeta[0] >= polygonal_vertex_arclength[v - 1].first) &&
5287  (zeta[0] <= polygonal_vertex_arclength[v].first))
5288  {
5289  // Increment in intrinsic coordinate along geom object
5290  double delta_zeta =
5291  (polygonal_vertex_arclength[v].second -
5292  polygonal_vertex_arclength[v - 1].second);
5293  // Increment in arclength along segment
5294  double delta_polyarc =
5295  (polygonal_vertex_arclength[v].first -
5296  polygonal_vertex_arclength[v - 1].first);
5297 
5298  // Mapped arclength coordinate
5299  double zeta_new =
5300  polygonal_vertex_arclength[v - 1].second +
5301  delta_zeta *
5302  (zeta[0] - polygonal_vertex_arclength[v - 1].first) /
5303  delta_polyarc;
5304  zeta[0] = zeta_new;
5305 
5306  // Success!
5307  found = true;
5308 
5309  // Bail out
5310  break;
5311  }
5312  } // for (v < n_vertex)
5313 
5314  // If we still haven't found it's probably the last point along
5315  if (!found)
5316  {
5317 #ifdef PARANOID
5318  double diff = std::fabs(
5319  zeta[0] - polygonal_vertex_arclength[n_vertex - 1].first);
5320  if (diff >
5321  ToleranceForVertexMismatchInPolygons::Tolerable_error)
5322  {
5323  std::ostringstream error_stream;
5324  error_stream
5325  << "Wasn't able to locate the polygonal arclength exactly\n"
5326  << "during re-setup of boundary coordinates and have\n"
5327  << "assumed that we're dealing with the final point along\n"
5328  << "the curvilinear segment and encountered some roundoff\n"
5329  << "However,the difference in the polygonal zeta "
5330  "coordinates\n"
5331  << "between zeta[0] " << zeta[0]
5332  << " and the originallly stored value "
5333  << polygonal_vertex_arclength[n_vertex - 1].first << "\n"
5334  << "is " << diff
5335  << " which exceeds the threshold specified\n"
5336  << "in the publically modifiable variable\n"
5337  << "ToleranceForVertexMismatchInPolygons::Tolerable_error\n"
5338  << "whose current value is: "
5339  << ToleranceForVertexMismatchInPolygons::Tolerable_error
5340  << "\nPlease check your mesh carefully and increase the\n"
5341  << "threshold if you're sure this is appropriate\n";
5342  throw OomphLibError(error_stream.str(),
5343  OOMPH_CURRENT_FUNCTION,
5344  OOMPH_EXCEPTION_LOCATION);
5345  }
5346 #endif
5347  zeta[0] = polygonal_vertex_arclength[n_vertex - 1].second;
5348  }
5349  }
5350 
5351  // Assign updated coordinate
5352  nod_pt->set_coordinates_on_boundary(b, zeta);
5353  }
5354  } // if (geom_object_pt != 0)
5355  else
5356  {
5357  // Establish the initial zeta value for this segment
5358  double z_initial = initial_segment_arclength[is];
5359 
5360  // Only use end points of the whole boundary and pick the
5361  // bottom left node
5362 
5363  // Set the bottom left coordinate as the first coordinates
5364  Vector<double> bottom_left_coordinate(2);
5365  bottom_left_coordinate = first_coordinate;
5366 
5367  // ... do the same with the zeta value
5368  Vector<double> bottom_left_zeta_coordinate(1);
5369  bottom_left_zeta_coordinate = first_node_zeta_coordinate;
5370 
5371  // Is the last y-coordinate smaller than y-coordinate of the
5372  // current bottom left coordinate
5373  if (last_coordinate[1] < bottom_left_coordinate[1])
5374  {
5375  // Re-set the bottom-left coordinate as the last coordinate
5376  bottom_left_coordinate = last_coordinate;
5377 
5378  // ... do the same with the zeta value
5379  bottom_left_zeta_coordinate = last_node_zeta_coordinate;
5380  }
5381  // The y-coordinates are the same
5382  else if (last_coordinate[1] == bottom_left_coordinate[1])
5383  {
5384  // Then check for the x-coordinate, which is the most-left
5385  if (last_coordinate[0] < bottom_left_coordinate[0])
5386  {
5387  // Re-set the bottom-left coordinate as the last
5388  // coordinate
5389  bottom_left_coordinate = last_coordinate;
5390 
5391  // ... do the same with the zeta value
5392  bottom_left_zeta_coordinate = last_node_zeta_coordinate;
5393  }
5394  } // else (The y-coordinates are the same)
5395 
5396  Vector<double> zeta(1, 0.0);
5397 
5398  // Now adjust boundary coordinate so that the bottom left
5399  // node has a boundary coordinate of 0.0 and that zeta
5400  // increases away from that point
5401  zeta = bottom_left_zeta_coordinate;
5402  const double zeta_ref = zeta[0];
5403 
5404  // Also get the maximum zeta value
5405  double zeta_max = 0.0;
5406  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
5407  it != all_nodes_pt.end();
5408  it++)
5409  {
5410  Node* nod_pt = (*it);
5411  nod_pt->get_coordinates_on_boundary(b, zeta);
5412 
5413 #ifdef OOMPH_HAS_MPI
5414 
5415  // If the mesh is distributed and the initial and final
5416  // zeta values for the segment have been assigned then
5417  // check if the segment is inverted, we need to consider
5418  // that to correctly assgin the zeta values for the segment
5419  if (this->is_mesh_distributed() &&
5420  Assigned_segments_initial_zeta_values[b])
5421  {
5422  // Is the segment inverted, if that is the case then
5423  // invert the zeta values
5424  if (boundary_segment_inverted(b)[is])
5425  {
5426  zeta[0] = segment_arclength[is] - zeta[0];
5427  } // if (boundary_segment_inverted(b)[is])
5428  }
5429 
5430 #endif // #ifdef OOMPH_HAS_MPI
5431 
5432  // Set the zeta value
5433  zeta[0] += z_initial;
5434  // Adjust the value based on the bottom-left condition
5435  zeta[0] -= zeta_ref;
5436  // If direction is reversed, then take absolute value
5437  if (zeta[0] < 0.0)
5438  {
5439  zeta[0] = std::fabs(zeta[0]);
5440  }
5441  // Get the max zeta value (we will use it to scale the
5442  // values to [0,1])
5443  if (zeta[0] > zeta_max)
5444  {
5445  zeta_max = zeta[0];
5446  }
5447  nod_pt->set_coordinates_on_boundary(b, zeta);
5448  } // Loop through the nodes in the segment (boundary)
5449 
5450 #ifdef OOMPH_HAS_MPI
5451  // After assigning boundary coordinates, BUT before scaling,
5452  // copy the initial and final segment arclengths so that we
5453  // know if the values go in increasing or decreasing
5454  // order. This will be used to identify the correct direction
5455  // of the segments in the new meshes created in the
5456  // adaptation method.
5457  if (this->is_mesh_distributed() &&
5458  Assigned_segments_initial_zeta_values[b])
5459  {
5460  // Get the first face element
5461  FiniteElement* first_seg_ele_pt = segment_sorted_ele_pt[is].front();
5462 
5463 #ifdef PARANOID
5464  // Check if the face element is nonhalo, it shouldn't, but
5465  // better check
5466  if (first_seg_ele_pt->is_halo())
5467  {
5468  std::ostringstream error_message;
5469  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
5470  "setup_boundary_coordinates()";
5471  error_message << "The first face element in the (" << is
5472  << ")-th segment"
5473  << " is halo\n";
5474  throw OomphLibError(
5475  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
5476  } // if (first_seg_ele_pt->is_halo())
5477 #endif
5478 
5479  // Number of nodes
5480  const unsigned nnod = first_seg_ele_pt->nnode();
5481 
5482  // Get the first node of the current segment
5483  Node* first_seg_node_pt = first_seg_ele_pt->node_pt(0);
5484  if (is_inverted[first_seg_ele_pt])
5485  {
5486  first_seg_node_pt = first_seg_ele_pt->node_pt(nnod - 1);
5487  }
5488 
5489  // Get the last face element
5490  FiniteElement* last_seg_ele_pt = segment_sorted_ele_pt[is].back();
5491 
5492 #ifdef PARANOID
5493  // Check if the face element is nonhalo, it shouldn't, but
5494  // better check
5495  if (last_seg_ele_pt->is_halo())
5496  {
5497  std::ostringstream error_message;
5498  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
5499  "setup_boundary_coordinates()";
5500  error_message << "The last face element in the (" << is
5501  << ")-th segment"
5502  << " is halo\n";
5503  throw OomphLibError(
5504  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
5505  } // if (last_seg_ele_pt->is_halo())
5506 #endif
5507 
5508  // Get the last node of the current segment
5509  Node* last_seg_node_pt = last_seg_ele_pt->node_pt(nnod - 1);
5510  if (is_inverted[last_seg_ele_pt])
5511  {
5512  last_seg_node_pt = last_seg_ele_pt->node_pt(0);
5513  }
5514 
5515  // Now get the first and last node boundary coordinate values
5516  Vector<double> first_seg_arclen(1);
5517  Vector<double> last_seg_arclen(1);
5518 
5519  first_seg_node_pt->get_coordinates_on_boundary(b, first_seg_arclen);
5520  last_seg_node_pt->get_coordinates_on_boundary(b, last_seg_arclen);
5521 
5522  // Update the initial and final segments arclength
5523  boundary_segment_initial_arclength(b)[is] = first_seg_arclen[0];
5524  boundary_segment_final_arclength(b)[is] = last_seg_arclen[0];
5525 
5526  // Update the initial and final coordinates
5527  Vector<double> updated_segment_initial_coord(2);
5528  Vector<double> updated_segment_final_coord(2);
5529  for (unsigned k = 0; k < 2; k++)
5530  {
5531  updated_segment_initial_coord[k] = first_seg_node_pt->x(k);
5532  updated_segment_final_coord[k] = last_seg_node_pt->x(k);
5533  }
5534 
5535  // Set the updated initial coordinate
5536  boundary_segment_initial_coordinate(b)[is] =
5537  updated_segment_initial_coord;
5538 
5539  // Set the updated final coordinate
5540  boundary_segment_final_coordinate(b)[is] =
5541  updated_segment_final_coord;
5542 
5543  } // if (this->is_mesh_distributed() &&
5544  // Assigned_segments_initial_zeta_values[b])
5545 #endif // #ifdef OOMPH_HAS_MPI
5546 
5547  // The max. value will be incorrect if we are working with
5548  // distributed meshes where the current boundary has been
5549  // split by the distribution process. Here correct the
5550  // maximum value
5551  if (zeta_max < boundary_arclength)
5552  {
5553  zeta_max = boundary_arclength;
5554  }
5555 
5556  // Scale all surface coordinates so that all the points be on
5557  // the range [0, 1]
5558  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
5559  it != all_nodes_pt.end();
5560  it++)
5561  {
5562  // Get the node
5563  Node* nod_pt = (*it);
5564 
5565  // Get the boundary coordinate
5566  nod_pt->get_coordinates_on_boundary(b, zeta);
5567 
5568  // Scale the value of the current node
5569  zeta[0] /= zeta_max;
5570 
5571  // Set the new scaled value
5572  nod_pt->set_coordinates_on_boundary(b, zeta);
5573  } // Loop over the nodes
5574 
5575  } // else if (geom_object_pt != 0)
5576 
5577  } // for (is < nsegments)
5578 
5579  // =================================================================
5580  // END: Scale the boundary coordinates based on whether the
5581  // boundary has an associated GeomObj or not
5582  // =================================================================
5583 
5584  // Cleanup
5585  for (unsigned e = 0; e < nel; e++)
5586  {
5587  delete face_el_pt[e];
5588  face_el_pt[e] = 0;
5589  }
5590 
5591  } // if (nel > 0), from the beginning of the method
5592 
5593  // Indicate that boundary coordinate has been set up
5594  Boundary_coordinate_exists[b] = true;
5595  }
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::is_mesh_distributed
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
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::UnstructuredTwoDMeshGeometryBase::Suppress_warning_about_regions_and_boundaries
static bool Suppress_warning_about_regions_and_boundaries
Public static flag to suppress warning; defaults to false.
Definition: unstructured_two_d_mesh_geometry_base.h:1741
oomph::UnstructuredTwoDMeshGeometryBase::boundary_element_in_region_pt
FiniteElement * boundary_element_in_region_pt(const unsigned &b, const unsigned &r, const unsigned &e) const
Return pointer to the e-th element adjacent to boundary b in region r.
Definition: unstructured_two_d_mesh_geometry_base.h:1893
oomph::UnstructuredTwoDMeshGeometryBase::boundary_geom_object_pt
std::map< unsigned, GeomObject * > & boundary_geom_object_pt()
Return direct access to the geometric object storage.
Definition: unstructured_two_d_mesh_geometry_base.h:1842
oomph::UnstructuredTwoDMeshGeometryBase::face_index_at_boundary_in_region
int face_index_at_boundary_in_region(const unsigned &b, const unsigned &r, const unsigned &e) const
Return face index of the e-th element adjacent to boundary b in region r.
Definition: unstructured_two_d_mesh_geometry_base.h:1911
oomph::UnstructuredTwoDMeshGeometryBase::nboundary_element_in_region
unsigned nboundary_element_in_region(const unsigned &b, const unsigned &r) const
Return the number of elements adjacent to boundary b in region r.
Definition: unstructured_two_d_mesh_geometry_base.h:1874
oomph::UnstructuredTwoDMeshGeometryBase::nregion
unsigned nregion()
Return the number of regions specified by attributes.
Definition: unstructured_two_d_mesh_geometry_base.h:1757
oomph::UnstructuredTwoDMeshGeometryBase::Polygonal_vertex_arclength_info
std::map< unsigned, Vector< std::pair< double, double > > > Polygonal_vertex_arclength_info
Definition: unstructured_two_d_mesh_geometry_base.h:2630
oomph::UnstructuredTwoDMeshGeometryBase::boundary_coordinate_limits
std::map< unsigned, Vector< double > > & boundary_coordinate_limits()
Definition: unstructured_two_d_mesh_geometry_base.h:1849
k
char char char int int * k
Definition: level2_impl.h:374
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
MergeRestartFiles.found
bool found
Definition: MergeRestartFiles.py:24
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
calibrate.error
int error
Definition: calibrate.py:297
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References b, oomph::Mesh::Boundary_coordinate_exists, boundary_coordinate_limits(), boundary_element_in_region_pt(), oomph::Mesh::boundary_element_pt(), boundary_geom_object_pt(), e(), calibrate::error, boost::multiprecision::fabs(), oomph::Mesh::face_index_at_boundary(), face_index_at_boundary_in_region(), MergeRestartFiles::found, oomph::Node::get_coordinates_on_boundary(), i, oomph::Mesh::is_mesh_distributed(), j, k, max, min, oomph::Mesh::nboundary_element(), nboundary_element_in_region(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), nregion(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Polygonal_vertex_arclength_info, oomph::GeomObject::position(), Region_attribute, oomph::Node::set_coordinates_on_boundary(), size, sqrt(), oomph::Global_string_for_annotation::string(), Suppress_warning_about_regions_and_boundaries, oomph::ToleranceForVertexMismatchInPolygons::Tolerable_error, v, oomph::Node::x(), and Eigen::zeta().

◆ snap_nodes_onto_geometric_objects()

void oomph::UnstructuredTwoDMeshGeometryBase::snap_nodes_onto_geometric_objects ( )
protected

Snap the boundary nodes onto any curvilinear boundaries defined by geometric objects

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

4021  {
4022  // Loop over all boundaries
4023  unsigned n_bound = this->nboundary();
4024  for (unsigned b = 0; b < n_bound; b++)
4025  {
4026  // Find the geometric object
4027  GeomObject* const geom_object_pt = this->boundary_geom_object_pt(b);
4028 
4029  // If there is one
4030  if (geom_object_pt != 0)
4031  {
4032  Vector<double> b_coord(1);
4033  Vector<double> new_x(2);
4034  const unsigned n_boundary_node = this->nboundary_node(b);
4035  for (unsigned n = 0; n < n_boundary_node; ++n)
4036  {
4037  // Get the boundary node and coordinates
4038  Node* const nod_pt = this->boundary_node_pt(b, n);
4039  nod_pt->get_coordinates_on_boundary(b, b_coord);
4040 
4041  // Get the position and time history according to the underlying
4042  // geometric object, assuming that it has the same timestepper
4043  // as the nodes....
4044  unsigned n_tvalues =
4045  1 + nod_pt->position_time_stepper_pt()->nprev_values();
4046  for (unsigned t = 0; t < n_tvalues; ++t)
4047  {
4048  // Get the position according to the underlying geometric object
4049  geom_object_pt->position(t, b_coord, new_x);
4050 
4051  // Move the node
4052  for (unsigned i = 0; i < 2; i++)
4053  {
4054  nod_pt->x(t, i) = new_x[i];
4055  }
4056  }
4057  }
4058  }
4059  }
4060  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
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::nboundary
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
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
plotPSD.t
t
Definition: plotPSD.py:36

References b, boundary_geom_object_pt(), oomph::Mesh::boundary_node_pt(), oomph::Node::get_coordinates_on_boundary(), i, n, oomph::Mesh::nboundary(), oomph::Mesh::nboundary_node(), oomph::TimeStepper::nprev_values(), oomph::GeomObject::position(), oomph::Node::position_time_stepper_pt(), plotPSD::t, and oomph::Node::x().

Member Data Documentation

◆ Allow_automatic_creation_of_vertices_on_boundaries

bool oomph::UnstructuredTwoDMeshGeometryBase::Allow_automatic_creation_of_vertices_on_boundaries
protected

Flag to indicate whether the automatic creation of vertices along the boundaries by Triangle is allowed

Referenced by disable_automatic_creation_of_vertices_on_boundaries(), enable_automatic_creation_of_vertices_on_boundaries(), is_automatic_creation_of_vertices_on_boundaries_allowed(), oomph::QuadFromTriangleMesh< ELEMENT >::QuadFromTriangleMesh(), and oomph::TriangleMesh< ELEMENT >::TriangleMesh().

◆ Boundary_coordinate_limits

std::map<unsigned, Vector<double> > oomph::UnstructuredTwoDMeshGeometryBase::Boundary_coordinate_limits
protected

Storage for the limits of the boundary coordinates defined by the use of geometric objects. Only used for curvilinear boundaries.

Referenced by boundary_coordinate_limits().

◆ Boundary_curve_section_pt

std::map<unsigned, TriangleMeshCurveSection*> oomph::UnstructuredTwoDMeshGeometryBase::Boundary_curve_section_pt
protected

A map that stores the associated curve section of the specified boundary id

Referenced by boundary_polyline_pt().

◆ Boundary_geom_object_pt

std::map<unsigned, GeomObject*> oomph::UnstructuredTwoDMeshGeometryBase::Boundary_geom_object_pt
protected

Storage for the geometric objects associated with any boundaries.

Referenced by boundary_geom_object_pt().

◆ Boundary_region_element_pt

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

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

Referenced by boundary_element_in_region_pt(), and nboundary_element_in_region().

◆ Extra_holes_coordinates

Vector<Vector<double> > oomph::UnstructuredTwoDMeshGeometryBase::Extra_holes_coordinates
protected

Storage for extra coordinates for holes.

◆ Face_index_region_at_boundary

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

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

Referenced by face_index_at_boundary_in_region().

◆ Free_curve_section_pt

std::set<TriangleMeshCurveSection*> oomph::UnstructuredTwoDMeshGeometryBase::Free_curve_section_pt
protected

A set that contains the curve sections created by this object therefore it is necessary to free their associated memory

Referenced by oomph::QuadFromTriangleMesh< ELEMENT >::~QuadFromTriangleMesh(), and oomph::TriangleMesh< ELEMENT >::~TriangleMesh().

◆ Free_open_curve_pt

std::set<TriangleMeshOpenCurve*> oomph::UnstructuredTwoDMeshGeometryBase::Free_open_curve_pt
protected

A set that contains the open curves created by this object therefore it is necessary to free their associated memory

Referenced by oomph::QuadFromTriangleMesh< ELEMENT >::~QuadFromTriangleMesh(), and oomph::TriangleMesh< ELEMENT >::~TriangleMesh().

◆ Free_polygon_pt

std::set<TriangleMeshPolygon*> oomph::UnstructuredTwoDMeshGeometryBase::Free_polygon_pt
protected

A set that contains the polygons created by this object therefore it is necessary to free their associated memory

Referenced by oomph::QuadFromTriangleMesh< ELEMENT >::~QuadFromTriangleMesh(), and oomph::TriangleMesh< ELEMENT >::~TriangleMesh().

◆ Internal_open_curve_pt

Vector<TriangleMeshOpenCurve*> oomph::UnstructuredTwoDMeshGeometryBase::Internal_open_curve_pt
protected

Vector of open polylines that define internal curves.

◆ Internal_polygon_pt

Vector<TriangleMeshPolygon*> oomph::UnstructuredTwoDMeshGeometryBase::Internal_polygon_pt
protected

Vector of polygons that define internal polygons.

◆ Nodes_on_boundary_pt

std::map<unsigned, std::set<Node*> > oomph::UnstructuredTwoDMeshGeometryBase::Nodes_on_boundary_pt
protected

Stores a pointer to a set with all the nodes related with a boundary

Referenced by nodes_on_boundary_pt().

◆ Outer_boundary_pt

Vector<TriangleMeshPolygon*> oomph::UnstructuredTwoDMeshGeometryBase::Outer_boundary_pt
protected

Polygon that defines outer boundaries.

◆ Polygonal_vertex_arclength_info

std::map<unsigned, Vector<std::pair<double, double> > > oomph::UnstructuredTwoDMeshGeometryBase::Polygonal_vertex_arclength_info
protected

Storage for pairs of doubles representing: .first: the arclength along the polygonal representation of the curviline .second: the corresponding intrinsic coordinate on the associated geometric object at which the vertices on the specified boundary are located. Only used for boundaries represented by geom objects.

Referenced by setup_boundary_coordinates().

◆ Region_attribute

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

Vector of attributes associated with the elements in each region.

Referenced by nregion_attribute(), region_attribute(), and setup_boundary_coordinates().

◆ Region_element_pt

std::map<unsigned, Vector<FiniteElement*> > oomph::UnstructuredTwoDMeshGeometryBase::Region_element_pt
protected

Vector of elements in each region differentiated by attribute (the key of the map is the attribute)

Referenced by nregion(), nregion_element(), and region_element_pt().

◆ Regions_coordinates

std::map<unsigned, Vector<double> > oomph::UnstructuredTwoDMeshGeometryBase::Regions_coordinates
protected

Storage for extra coordinates for regions. The key on the map is the region id

◆ Suppress_warning_about_regions_and_boundaries

bool oomph::UnstructuredTwoDMeshGeometryBase::Suppress_warning_about_regions_and_boundaries
static

Public static flag to suppress warning; defaults to false.

Referenced by setup_boundary_coordinates().

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