triangle_mesh.template.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_TRIANGLE_MESH_TEMPLATE_CC
27 #define OOMPH_TRIANGLE_MESH_TEMPLATE_CC
28 
29 #include <iostream>
30 
31 #include "triangle_mesh.template.h"
32 #include "../generic/map_matrix.h"
33 #include "../generic/multi_domain.h"
34 #include "../generic/projection.h"
35 #include "../generic/face_element_as_geometric_object.h"
36 
37 namespace oomph
38 {
39  //======================================================================
40  /// Build with the help of the scaffold mesh coming
41  /// from the triangle mesh generator Triangle.
42  //======================================================================
43  template<class ELEMENT>
45  const bool& use_attributes)
46  {
47  // Mesh can only be built with 2D Telements.
48  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
49 
50  // Create space for elements
51  unsigned nelem = Tmp_mesh_pt->nelement();
52  Element_pt.resize(nelem);
53 
54  // Create space for nodes
55  unsigned nnode_scaffold = Tmp_mesh_pt->nnode();
56 
57  // Create a map storing the node_id of the mesh used to update the
58  // node position in the update_triangulateio function
59  std::map<Node*, unsigned> old_global_number;
60 
61  // Store the TriangulateIO node id
62  for (unsigned inod = 0; inod < nnode_scaffold; inod++)
63  {
64  Node* old_node_pt = Tmp_mesh_pt->node_pt(inod);
65  old_global_number[old_node_pt] = inod;
66  }
67 
68  // Initialize the old node id vector
69  Oomph_vertex_nodes_id.resize(nnode_scaffold);
70 
71  // Create space for nodes
72  Node_pt.resize(nnode_scaffold, 0);
73 
74  // Set number of boundaries
75  unsigned nbound = Tmp_mesh_pt->nboundary();
76 
77  // Resize the boundary information
78  set_nboundary(nbound);
79  Boundary_element_pt.resize(nbound);
80  Face_index_at_boundary.resize(nbound);
81 
82  // If we have different regions, then resize the region
83  // information
84  if (use_attributes)
85  {
86  Boundary_region_element_pt.resize(nbound);
87  Face_index_region_at_boundary.resize(nbound);
88  }
89 
90  // Loop over elements in scaffold mesh, visit their nodes
91  for (unsigned e = 0; e < nelem; e++)
92  {
93  Element_pt[e] = new ELEMENT;
94  }
95 
96  // Number of nodes per element from the scaffold mesh
97  unsigned nnod_el = Tmp_mesh_pt->finite_element_pt(0)->nnode();
98 
99  // Setup map to check the (pseudo-)global node number
100  // Nodes whose number is zero haven't been copied across
101  // into the mesh yet.
102  std::map<Node*, unsigned> global_number;
103  unsigned global_count = 0;
104 
105  // Map of Element attribute pairs
106  std::map<double, Vector<FiniteElement*>> element_attribute_map;
107 
108  // If we're using attributes
109  if (use_attributes)
110  {
111  // If we're using attributes then we need attribute 0 which will
112  // be associated with region 0
113  element_attribute_map[0].resize(0);
114  }
115 
116  // Loop over elements in scaffold mesh, visit their nodes
117  for (unsigned e = 0; e < nelem; e++)
118  {
119  // Loop over all nodes in element
120  for (unsigned j = 0; j < nnod_el; j++)
121  {
122  // Pointer to node in the scaffold mesh
123  Node* scaffold_node_pt = Tmp_mesh_pt->finite_element_pt(e)->node_pt(j);
124 
125  // Get the (pseudo-)global node number in scaffold mesh
126  // (It's zero [=default] if not visited this one yet)
127  unsigned j_global = global_number[scaffold_node_pt];
128 
129  // Haven't done this one yet
130  if (j_global == 0)
131  {
132  // Find and store the node_id in the old nodes map
133  Oomph_vertex_nodes_id[global_count] =
134  old_global_number[scaffold_node_pt];
135 
136  // Get pointer to set of mesh boundaries that this
137  // scaffold node occupies; NULL if the node is not on any boundary
138  std::set<unsigned>* boundaries_pt;
139  scaffold_node_pt->get_boundaries_pt(boundaries_pt);
140 
141  // Storage for the new node
142  Node* new_node_pt = 0;
143 
144  // Is it on boundaries
145  if (boundaries_pt != 0)
146  {
147  // Create new boundary node
148  new_node_pt =
149  finite_element_pt(e)->construct_boundary_node(j, time_stepper_pt);
150 
151  // Add to boundaries
152  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
153  it != boundaries_pt->end();
154  ++it)
155  {
156  add_boundary_node(*it, new_node_pt);
157  }
158  }
159  // Build normal node
160  else
161  {
162  // Create new normal node
163  new_node_pt =
164  finite_element_pt(e)->construct_node(j, time_stepper_pt);
165  }
166 
167  // Give it a number (not necessarily the global node
168  // number in the scaffold mesh -- we just need something
169  // to keep track...)
170  global_count++;
171  global_number[scaffold_node_pt] = global_count;
172 
173  // Copy new node, created using the NEW element's construct_node
174  // function into global storage, using the same global
175  // node number that we've just associated with the
176  // corresponding node in the scaffold mesh
177  Node_pt[global_count - 1] = new_node_pt;
178 
179  // Assign coordinates
180  for (unsigned i = 0; i < finite_element_pt(e)->dim(); i++)
181  {
182  new_node_pt->x(i) = scaffold_node_pt->x(i);
183  }
184  }
185  // This one has already been done: Copy accross
186  else
187  {
188  finite_element_pt(e)->node_pt(j) = Node_pt[j_global - 1];
189  }
190  }
191 
192  // If we're using attributes
193  if (use_attributes)
194  {
195  element_attribute_map[Tmp_mesh_pt->element_attribute(e)].push_back(
196  finite_element_pt(e));
197  }
198  }
199 
200  // Now let's construct lists
201  // Find the number of attributes
202  if (use_attributes)
203  {
204  unsigned n_attribute = element_attribute_map.size();
205 
206  // There are n_attribute different regions
207  this->Region_attribute.resize(n_attribute);
208 
209  // Copy the vectors in the map over to our internal storage
210  unsigned count = 0;
211  for (std::map<double, Vector<FiniteElement*>>::iterator it =
212  element_attribute_map.begin();
213  it != element_attribute_map.end();
214  ++it)
215  {
216  this->Region_attribute[count] = it->first;
217  Region_element_pt[static_cast<unsigned>(Region_attribute[count])] =
218  it->second;
219  ++count;
220  }
221  }
222 
223  // At this point we've created all the elements and
224  // created their vertex nodes. Now we need to create
225  // the additional (midside and internal) nodes!
226 
227  unsigned boundary_id = 0;
228 
229  // Get number of nodes along element edge and dimension of element (2)
230  // from first element
231  unsigned n_node_1d = finite_element_pt(0)->nnode_1d();
232  unsigned dim = finite_element_pt(0)->dim();
233 
234  // Storage for the local coordinate of the new node
235  Vector<double> s(dim);
236 
237  // Get number of nodes in the element from first element
238  unsigned n_node = finite_element_pt(0)->nnode();
239 
240  // Storage for each global edge of the mesh
241  unsigned n_global_edge = Tmp_mesh_pt->nglobal_edge();
242  Vector<Vector<Node*>> nodes_on_global_edge(n_global_edge);
243 
244  // Loop over elements
245  for (unsigned e = 0; e < nelem; e++)
246  {
247  // Cache pointers to the elements
248  FiniteElement* const elem_pt = finite_element_pt(e);
249  FiniteElement* const tmp_elem_pt = Tmp_mesh_pt->finite_element_pt(e);
250 
251  // The number of edge nodes is 3*(nnode_1d-1)
252  unsigned n_edge_node = 3 * (n_node_1d - 1);
253 
254  // If there are any more nodes, these are internal and can be
255  // constructed and added directly to the mesh
256  for (unsigned n = n_edge_node; n < n_node; ++n)
257  {
258  // Create new node (it can never be a boundary node)
259  Node* new_node_pt = elem_pt->construct_node(n, time_stepper_pt);
260 
261  // What are the node's local coordinates?
262  elem_pt->local_coordinate_of_node(n, s);
263 
264  // Find the coordinates of the new node from the existing
265  // and fully-functional element in the scaffold mesh
266  for (unsigned i = 0; i < dim; i++)
267  {
268  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s, i);
269  }
270 
271  // Add the node to the mesh's global look-up scheme
272  Node_pt.push_back(new_node_pt);
273  }
274 
275  // Now loop over the mid-side edge nodes
276  // Start from node number 3
277  unsigned n = 3;
278 
279  // Loop over edges
280  for (unsigned j = 0; j < 3; j++)
281  {
282  // Find the boundary id of the edge
283  boundary_id = Tmp_mesh_pt->edge_boundary(e, j);
284 
285  // Find the global edge index
286  unsigned edge_index = Tmp_mesh_pt->edge_index(e, j);
287 
288  // If the nodes on the edge have not been allocated, construct them
289  if (nodes_on_global_edge[edge_index].size() == 0)
290  {
291  // Loop over the nodes on the edge excluding the ends
292  for (unsigned j2 = 0; j2 < n_node_1d - 2; ++j2)
293  {
294  // Storage for the new node
295  Node* new_node_pt = 0;
296 
297  // If the edge is on a boundary, construct a boundary node
298  if (boundary_id > 0)
299  {
300  new_node_pt =
301  elem_pt->construct_boundary_node(n, time_stepper_pt);
302  // Add it to the boundary
303  this->add_boundary_node(boundary_id - 1, new_node_pt);
304  }
305  // Otherwise construct a normal node
306  else
307  {
308  new_node_pt = elem_pt->construct_node(n, time_stepper_pt);
309  }
310 
311  // What are the node's local coordinates?
312  elem_pt->local_coordinate_of_node(n, s);
313 
314  // Find the coordinates of the new node from the existing
315  // and fully-functional element in the scaffold mesh
316  for (unsigned i = 0; i < dim; i++)
317  {
318  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s, i);
319  }
320 
321  // Add to the global node list
322  Node_pt.push_back(new_node_pt);
323 
324  // Add to the edge index
325  nodes_on_global_edge[edge_index].push_back(new_node_pt);
326  // Increment the node number
327  ++n;
328  }
329  }
330  // Otherwise just set the pointers
331  // using the fact that the next time the edge is visited
332  // the nodes must be arranged in the other order because all
333  // triangles have the same orientation
334  else
335  {
336  // Loop over the nodes on the edge excluding the ends
337  for (unsigned j2 = 0; j2 < n_node_1d - 2; ++j2)
338  {
339  // Set the local node from the edge but indexed the other
340  // way around
341  elem_pt->node_pt(n) =
342  nodes_on_global_edge[edge_index][n_node_1d - 3 - j2];
343  ++n;
344  }
345  }
346 
347  // Set the elements adjacent to the boundary from the
348  // boundary id information
349  if (boundary_id > 0)
350  {
351  Boundary_element_pt[boundary_id - 1].push_back(elem_pt);
352  // Need to put a shift in here because of an inconsistent naming
353  // convention between triangle and face elements
354  Face_index_at_boundary[boundary_id - 1].push_back((j + 2) % 3);
355 
356  // If using regions set up the boundary information
357  if (use_attributes)
358  {
359  unsigned tmp_region =
360  static_cast<unsigned>(Tmp_mesh_pt->element_attribute(e));
361  // Element adjacent to boundary
362  Boundary_region_element_pt[boundary_id - 1][tmp_region].push_back(
363  elem_pt);
364  // Need to put a shift in here because of an inconsistent naming
365  // convention between triangle and face elements
366  Face_index_region_at_boundary[boundary_id - 1][tmp_region]
367  .push_back((j + 2) % 3);
368  }
369  }
370 
371  } // end of loop over edges
372  } // end of loop over elements
373 
374 
375  // Lookup scheme has now been setup
376  Lookup_for_elements_next_boundary_is_setup = true;
377  }
378 
379 #ifdef OOMPH_HAS_MPI
380 
381  //======================================================================
382  /// Identify the segments from the old mesh (original mesh)
383  /// in the new mesh (this) and assign initial and final boundary
384  /// coordinates for the segments that create the boundary
385  //======================================================================
386  template<class ELEMENT>
389  const unsigned& b, TriangleMesh<ELEMENT>* original_mesh_pt)
390  {
391  // ------------------------------------------------------------------
392  // First: Get the face elements associated with the current boundary
393  // (nonhalo elements only)
394  // ------------------------------------------------------------------
395  // Temporary storage for face elements
396  Vector<FiniteElement*> face_el_pt;
397 
398  // Temporary storage for number of elements adjacent to the boundary
399  unsigned nele = 0;
400 
401  // Temporary storage for elements adjacent to the boundary that have
402  // a common edge (related with internal boundaries)
403  unsigned n_repeated_ele = 0;
404 
405  const unsigned n_regions = this->nregion();
406 
407  // map to associate the face element to the bulk element, necessary
408  // to attach halo face elements at both sides of each found segment
409  std::map<FiniteElement*, FiniteElement*> face_to_bulk_element_pt;
410 
411  // Temporary storage for already done nodes
412  Vector<std::pair<Node*, Node*>> done_nodes_pt;
413 
414  // If there is more than one region then only use boundary
415  // coordinates from the bulk side (region 0)
416  if (n_regions > 1)
417  {
418  for (unsigned rr = 0; rr < n_regions; rr++)
419  {
420  const unsigned region_id =
421  static_cast<unsigned>(this->Region_attribute[rr]);
422 
423  // Loop over all elements on boundaries in region i_r
424  const unsigned nel_in_region =
425  this->nboundary_element_in_region(b, region_id);
426 
427  unsigned nel_repetead_in_region = 0;
428 
429  // Only bother to do anything else, if there are elements
430  // associated with the boundary and the current region
431  if (nel_in_region > 0)
432  {
433  // Flag that activates when a repeated face element is found,
434  // possibly because we are dealing with an internal boundary
435  bool repeated = false;
436 
437  // Loop over the bulk elements adjacent to boundary b
438  for (unsigned e = 0; e < nel_in_region; e++)
439  {
440  // Get pointer to the bulk element that is adjacent to boundary b
441  FiniteElement* bulk_elem_pt =
442  this->boundary_element_in_region_pt(b, region_id, e);
443 
444 #ifdef OOMPH_HAS_MPI
445  // In a distributed mesh only work with nonhalo elements
446  if (this->is_mesh_distributed() && bulk_elem_pt->is_halo())
447  {
448  // Increase the number of repeated elements
449  n_repeated_ele++;
450  // Go for the next element
451  continue;
452  }
453 #endif
454 
455  // Find the index of the face of element e along boundary b
456  int face_index =
457  this->face_index_at_boundary_in_region(b, region_id, e);
458 
459  // Before adding the new element we need to be sure that
460  // the edge that this element represent has not been
461  // already added
462  FiniteElement* tmp_ele_pt =
463  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
464 
465  const unsigned n_nodes = tmp_ele_pt->nnode();
466 
467  std::pair<Node*, Node*> tmp_pair = std::make_pair(
468  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
469 
470  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
471  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
472 
473  // Search for repeated nodes
474  const unsigned n_done_nodes = done_nodes_pt.size();
475  for (unsigned l = 0; l < n_done_nodes; l++)
476  {
477  if (tmp_pair == done_nodes_pt[l] ||
478  tmp_pair_inverse == done_nodes_pt[l])
479  {
480  nel_repetead_in_region++;
481  repeated = true;
482  break;
483  }
484  }
485 
486  // Create new face element
487  if (!repeated)
488  {
489  // Add the pair of nodes (edge) to the node dones
490  done_nodes_pt.push_back(tmp_pair);
491  // Create the map to know if the element is halo
492  face_el_pt.push_back(tmp_ele_pt);
493  // Add the element to the face elements
494  face_to_bulk_element_pt[tmp_ele_pt] = bulk_elem_pt;
495  }
496  else
497  {
498  // Clean up
499  delete tmp_ele_pt;
500  tmp_ele_pt = 0;
501  }
502 
503  // Re-start
504  repeated = false;
505 
506  } // for (e < nel_in_region)
507 
508  nele += nel_in_region;
509 
510  n_repeated_ele += nel_repetead_in_region;
511 
512  } // if (nel_in_region > 0)
513  } // for (rr < n_regions)
514  } // if (n_regions > 1)
515  // Otherwise it's just the normal boundary functions
516  else
517  {
518  // Loop over all elements on boundaries
519  nele = this->nboundary_element(b);
520 
521  // Only bother to do anything else, if there are elements
522  if (nele > 0)
523  {
524  // Flag that activates when a repeated face element is found,
525  // possibly because we are dealing with an internal boundary
526  bool repeated = false;
527 
528  // Loop over the bulk elements adjacent to boundary b
529  for (unsigned e = 0; e < nele; e++)
530  {
531  // Get pointer to the bulk element that is adjacent to boundary b
532  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
533 
534 #ifdef OOMPH_HAS_MPI
535  // In a distributed mesh only work with nonhalo elements
536  if (this->is_mesh_distributed() && bulk_elem_pt->is_halo())
537  {
538  // Increase the number of repeated elements
539  n_repeated_ele++;
540  // Go for the next element
541  continue;
542  }
543 #endif
544 
545  // Find the index of the face of element e along boundary b
546  int face_index = this->face_index_at_boundary(b, e);
547 
548  // Before adding the new element we need to be sure that
549  // the edge that this element represents has not been
550  // already added (only applies for internal boundaries)
551  FiniteElement* tmp_ele_pt =
552  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
553 
554  const unsigned n_nodes = tmp_ele_pt->nnode();
555 
556  std::pair<Node*, Node*> tmp_pair = std::make_pair(
557  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
558 
559  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
560  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
561 
562  // Search for repeated nodes
563  const unsigned n_done_nodes = done_nodes_pt.size();
564  for (unsigned l = 0; l < n_done_nodes; l++)
565  {
566  if (tmp_pair == done_nodes_pt[l] ||
567  tmp_pair_inverse == done_nodes_pt[l])
568  {
569  // Increase the number of repeated elements
570  n_repeated_ele++;
571  // Mark the element as repeated
572  repeated = true;
573  break;
574  }
575  }
576 
577  // Create new face element
578  if (!repeated)
579  {
580  // Add the pair of nodes (edge) to the node dones
581  done_nodes_pt.push_back(tmp_pair);
582  // Add the element to the face elements
583  face_el_pt.push_back(tmp_ele_pt);
584  // Create the map to know if the element is halo
585  face_to_bulk_element_pt[tmp_ele_pt] = bulk_elem_pt;
586  }
587  else
588  {
589  // Free the repeated bulk element!!
590  delete tmp_ele_pt;
591  tmp_ele_pt = 0;
592  }
593 
594  // Re-start
595  repeated = false;
596 
597  } // for (e < nel)
598  } // if (nel > 0)
599 
600  } // else (n_regions > 1)
601 
602  // Do not consider the repeated elements
603  nele -= n_repeated_ele;
604 
605 #ifdef PARANOID
606  if (nele != face_el_pt.size())
607  {
608  std::ostringstream error_message;
609  error_message
610  << "The independent counting of face elements (" << nele << ") for "
611  << "boundary (" << b << ") is different\n"
612  << "from the real number of face elements in the container ("
613  << face_el_pt.size() << ")\n";
614  throw OomphLibError(error_message.str(),
615  "TriangleMesh::identify_boundary_segments_and_assign_"
616  "initial_zeta_values()",
617  OOMPH_EXCEPTION_LOCATION);
618  }
619 #endif
620 
621  // Continue even thought there are no elements, the processor needs
622  // to participate in the communications
623 
624  // ----------------------------------------------------------------
625  // Second: Sort the face elements, only consider nonhalo elements
626  // ----------------------------------------------------------------
627 
628  // A flag vector to mark those face elements that are considered as
629  // halo in the current processor
630  std::vector<bool> is_halo_face_element(nele, false);
631 
632  // Count the total number of non halo face elements
633  unsigned nnon_halo_face_elements = 0;
634 
635  // We will have halo face elements if the mesh is distributed
636  for (unsigned ie = 0; ie < nele; ie++)
637  {
638  // Get the face element
639  FiniteElement* face_ele_pt = face_el_pt[ie];
640  // Get the bulk element
641  FiniteElement* tmp_bulk_ele_pt = face_to_bulk_element_pt[face_ele_pt];
642  // Check if the bulk element is halo
643  if (!tmp_bulk_ele_pt->is_halo())
644  {
645  is_halo_face_element[ie] = false;
646  nnon_halo_face_elements++;
647  }
648  else
649  {
650  // Mark the face element as halo
651  is_halo_face_element[ie] = true;
652  }
653  } // for (ie < nele)
654 
655 #ifdef PARANOID
656  // Get the total number of halo face elements
657  const unsigned nhalo_face_element = nele - nnon_halo_face_elements;
658  if (nhalo_face_element > 0)
659  {
660  std::ostringstream error_message;
661  error_message
662  << "There should not be halo face elements since they were not "
663  << "considered when computing the face elements\n\n"
664  << "The number of found halo face elements is: " << nhalo_face_element
665  << "\n\n";
666  throw OomphLibError(error_message.str(),
667  "TriangleMesh::identify_boundary_segments_and_assign_"
668  "initial_zeta_values()",
669  OOMPH_EXCEPTION_LOCATION);
670  }
671 #endif
672 
673  // The vector of list to store the "segments" that compound the
674  // boundary (segments may appear only in a distributed mesh)
675  Vector<std::list<FiniteElement*>> segment_sorted_ele_pt;
676 
677  // Number of already sorted face elements (only nonhalo elements for
678  // a distributed mesh)
679  unsigned nsorted_face_elements = 0;
680 
681  // Keep track of who's done (this apply to nonhalo only, remember we
682  // are only working with nonhalo elements)
683  std::map<FiniteElement*, bool> done_el;
684 
685  // Keep track of which element is inverted (in distributed mesh the
686  // elements may be inverted with respect to the segment they belong)
687  std::map<FiniteElement*, bool> is_inverted;
688 
689  // Iterate until all possible segments have been created
690  while (nsorted_face_elements < nnon_halo_face_elements)
691  {
692  // The ordered list of face elements (in a distributed mesh a
693  // collection of contiguous face elements define a segment)
694  std::list<FiniteElement*> sorted_el_pt;
695  sorted_el_pt.clear();
696 
697 #ifdef PARANOID
698  // Select an initial element for the segment
699  bool found_initial_face_element = false;
700 #endif
701 
702  FiniteElement* ele_face_pt = 0;
703 
704  unsigned iface = 0;
705  for (iface = 0; iface < nele; iface++)
706  {
707  if (!is_halo_face_element[iface])
708  {
709  ele_face_pt = face_el_pt[iface];
710  // If not done then take it as initial face element
711  if (!done_el[ele_face_pt])
712  {
713 #ifdef PARANOID
714  found_initial_face_element = true;
715 #endif
716  nsorted_face_elements++;
717  iface++; // The next element number
718  sorted_el_pt.push_back(ele_face_pt);
719  // Mark as done
720  done_el[ele_face_pt] = true;
721  break;
722  }
723  }
724  } // for (iface < nele)
725 
726 #ifdef PARANOID
727  if (!found_initial_face_element)
728  {
729  std::ostringstream error_message;
730  error_message
731  << "Could not find an initial face element for the current segment\n";
732  throw OomphLibError(error_message.str(),
733  "TriangleMesh::identify_boundary_segments_and_"
734  "assign_initial_zeta_values()",
735  OOMPH_EXCEPTION_LOCATION);
736  }
737 #endif
738 
739  // Number of nodes
740  const unsigned nnod = ele_face_pt->nnode();
741 
742  // Left and right most nodes (the left and right nodes of the
743  // current face element)
744  Node* left_node_pt = ele_face_pt->node_pt(0);
745  Node* right_node_pt = ele_face_pt->node_pt(nnod - 1);
746 
747  // Continue iterating if a new face element has been added to the
748  // list
749  bool face_element_added = false;
750 
751  // While a new face element has been added to the set of sorted
752  // face elements then re-iterate
753  do
754  {
755  // Start from the next face element since we have already added
756  // the previous one as the initial face element (any previous
757  // face element had to be added on previous iterations)
758  for (unsigned iiface = iface; iiface < nele; iiface++)
759  {
760  // Re-start flag
761  face_element_added = false;
762 
763  // Get the candidate element
764  ele_face_pt = face_el_pt[iiface];
765 
766  // Check that the candidate element has not been done and is
767  // not a halo element
768  if (!(done_el[ele_face_pt] || is_halo_face_element[iiface]))
769  {
770  // Get the left and right nodes of the current element
771  Node* local_left_node_pt = ele_face_pt->node_pt(0);
772  Node* local_right_node_pt = ele_face_pt->node_pt(nnod - 1);
773  // New element fits at the left of segment and is not inverted
774  if (left_node_pt == local_right_node_pt)
775  {
776  left_node_pt = local_left_node_pt;
777  sorted_el_pt.push_front(ele_face_pt);
778  is_inverted[ele_face_pt] = false;
779  face_element_added = true;
780  }
781  // New element fits at the left of segment and is inverted
782  else if (left_node_pt == local_left_node_pt)
783  {
784  left_node_pt = local_right_node_pt;
785  sorted_el_pt.push_front(ele_face_pt);
786  is_inverted[ele_face_pt] = true;
787  face_element_added = true;
788  }
789  // New element fits on the right of segment and is not inverted
790  else if (right_node_pt == local_left_node_pt)
791  {
792  right_node_pt = local_right_node_pt;
793  sorted_el_pt.push_back(ele_face_pt);
794  is_inverted[ele_face_pt] = false;
795  face_element_added = true;
796  }
797  // New element fits on the right of segment and is inverted
798  else if (right_node_pt == local_right_node_pt)
799  {
800  right_node_pt = local_left_node_pt;
801  sorted_el_pt.push_back(ele_face_pt);
802  is_inverted[ele_face_pt] = true;
803  face_element_added = true;
804  }
805 
806  if (face_element_added)
807  {
808  done_el[ele_face_pt] = true;
809  nsorted_face_elements++;
810  break;
811  }
812 
813  } // if (!(done_el[ele_face_pt] || is_halo_face_element[iiface]))
814  } // for (iiface<nnon_halo_face_element)
815  } while (face_element_added &&
816  (nsorted_face_elements < nnon_halo_face_elements));
817 
818  // Store the created segment in the vector of segments
819  segment_sorted_ele_pt.push_back(sorted_el_pt);
820 
821  } // while(nsorted_face_elements < nnon_halo_face_elements);
822 
823  // The number of segments in this processor
824  const unsigned nsegments = segment_sorted_ele_pt.size();
825 
826  // ------------------------------------------------------------------
827  // Third: We have the face elements sorted (nonhalo only), now
828  // assign boundary coordinates to the nodes in the segments. This is
829  // the LOCAL boundary coordinate which is required if the zeta
830  // values need to be inverted
831  // ------------------------------------------------------------------
832  // Necessary in case boundaries with no geom object associated need
833  // to be inverted the zeta values (It is necessary to compute the
834  // arclength but also to store the nodes in a container (set))
835  // ------------------------------------------------------------------
836 
837  // Vector of sets that stores the nodes of each segment based on a
838  // lexicographically order starting from the bottom left node of
839  // each segment
840  Vector<std::set<Node*>> segment_all_nodes_pt;
841 
842  // The arclength of each segment in the current processor
843  Vector<double> segment_arclength(nsegments);
844 
845  // The number of vertices of each segment
846  Vector<unsigned> nvertices_per_segment(nsegments);
847 
848  // The initial zeta for the segment
849  Vector<double> initial_zeta_segment(nsegments);
850 
851  // The final zeta for the segment
852  Vector<double> final_zeta_segment(nsegments);
853 
854 #ifdef PARANOID
855  if (nnon_halo_face_elements > 0 && nsegments == 0)
856  {
857  std::ostringstream error_message;
858  error_message
859  << "The number of segments is zero, but the number of nonhalo\n"
860  << "elements is: (" << nnon_halo_face_elements << ")\n";
861  throw OomphLibError(error_message.str(),
862  "TriangleMesh::identify_boundary_segments_and_assign_"
863  "initial_zeta_values()",
864  OOMPH_EXCEPTION_LOCATION);
865  } // if (nnon_halo_face_elements > 0 && nsegments == 0)
866 #endif
867 
868  // Go through all the segments and compute the LOCAL boundary
869  // coordinates
870  for (unsigned is = 0; is < nsegments; is++)
871  {
872 #ifdef PARANOID
873  if (segment_sorted_ele_pt[is].size() == 0)
874  {
875  std::ostringstream error_message;
876  error_message << "The (" << is << ")-th segment has no elements\n";
877  throw OomphLibError(error_message.str(),
878  "TriangleMesh::identify_boundary_segments_and_"
879  "assign_initial_zeta_values()",
880  OOMPH_EXCEPTION_LOCATION);
881  } // if (segment_sorted_ele_pt[is].size() == 0)
882 #endif
883 
884  // Get access to the first element on the segment
885  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
886 
887  // Number of nodes
888  const unsigned nnod = first_ele_pt->nnode();
889 
890  // Get the first node of the current segment
891  Node* first_node_pt = first_ele_pt->node_pt(0);
892  if (is_inverted[first_ele_pt])
893  {
894  first_node_pt = first_ele_pt->node_pt(nnod - 1);
895  }
896 
897  // Get access to the last element on the segment
898  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
899 
900  // Get the last node of the current segment
901  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
902  if (is_inverted[last_ele_pt])
903  {
904  last_node_pt = last_ele_pt->node_pt(0);
905  }
906 
907  // Coordinates of left node
908  double x_left = first_node_pt->x(0);
909  double y_left = first_node_pt->x(1);
910 
911  // Initialise boundary coordinate (local boundary coordinate for
912  // boundaries with more than one segment)
913  Vector<double> zeta(1, 0.0);
914 
915  // If the boundary has an associated GeomObject then it is not
916  // necessary to compute the arclength, only read the values from
917  // the nodes at the edges
918  if (this->boundary_geom_object_pt(b) != 0)
919  {
920  first_node_pt->get_coordinates_on_boundary(b, zeta);
921  initial_zeta_segment[is] = zeta[0];
922  last_node_pt->get_coordinates_on_boundary(b, zeta);
923  final_zeta_segment[is] = zeta[0];
924  }
925 
926  // Lexicographically bottom left node
927  std::set<Node*> local_nodes_pt;
928  local_nodes_pt.insert(first_node_pt);
929 
930  // Now loop over nodes in order
931  for (std::list<FiniteElement*>::iterator it =
932  segment_sorted_ele_pt[is].begin();
933  it != segment_sorted_ele_pt[is].end();
934  it++)
935  {
936  // Get element
937  FiniteElement* el_pt = *it;
938 
939  // Start node and increment
940  unsigned k_nod = 1;
941  int nod_diff = 1;
942  if (is_inverted[el_pt])
943  {
944  k_nod = nnod - 2;
945  nod_diff = -1;
946  }
947 
948  // Loop over nodes
949  for (unsigned j = 1; j < nnod; j++)
950  {
951  Node* nod_pt = el_pt->node_pt(k_nod);
952  k_nod += nod_diff;
953 
954  // Coordinates of right node
955  double x_right = nod_pt->x(0);
956  double y_right = nod_pt->x(1);
957 
958  // Increment boundary coordinate (the arclength)
959  zeta[0] += sqrt((x_right - x_left) * (x_right - x_left) +
960  (y_right - y_left) * (y_right - y_left));
961 
962  // // When we have a GeomObject associated to the boundary we already
963  // // know the zeta values for the nodes, there is no need to compute
964  // // the arclength
965  // if (this->boundary_geom_object_pt(b)==0)
966  // {
967  // // Set boundary coordinate
968  // nod_pt->set_coordinates_on_boundary(b, zeta);
969  // }
970 
971  // Increment reference coordinate
972  x_left = x_right;
973  y_left = y_right;
974 
975  // Get lexicographically bottom left node but only
976  // use vertex nodes as candidates
977  local_nodes_pt.insert(nod_pt);
978  } // for (j < nnod)
979 
980  } // iterator over the elements in the segment
981 
982  // Store the arclength of the segment
983  segment_arclength[is] = zeta[0];
984 
985  // Store the number of vertices in the segment
986  nvertices_per_segment[is] = local_nodes_pt.size();
987 
988  // Add the nodes for the corresponding segment in the container
989  segment_all_nodes_pt.push_back(local_nodes_pt);
990 
991  } // for (is < nsegments)
992 
993  // Get the number of sets for nodes
994 #ifdef PARANOID
995  if (segment_all_nodes_pt.size() != nsegments)
996  {
997  std::ostringstream error_message;
998  error_message << "The number of segments (" << nsegments
999  << ") and the number of "
1000  << "sets of nodes (" << segment_all_nodes_pt.size()
1001  << ") representing\n"
1002  << "the\nsegments is different!!!\n\n";
1003  throw OomphLibError(
1004  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1005  }
1006 #endif
1007 
1008  // Store the initial arclength for each segment of boundary in the
1009  // current processor, initalise to zero in case we have a non
1010  // distributed boundary
1011  Vector<double> initial_segment_arclength(nsegments, 0.0);
1012 
1013  // Associated the index of the current segment to the segment index
1014  // in the original mesh (input mesh)
1015  Vector<unsigned> current_segment_to_original_segment_index(nsegments);
1016 
1017  // Each segment needs to know whether it has to be inverted or not
1018  // Store whether a segment needs to be inverted or not
1019  Vector<unsigned> segment_inverted(nsegments);
1020 
1021  // -----------------------------------------------------------------
1022  // Fourth: Identify the segments with the ones in the original mesh
1023  // (has sense only in the adaptation process)
1024  // -----------------------------------------------------------------
1025 
1026  // Now check if there are segments associated to this boundary
1027  if (nsegments > 0)
1028  {
1029 #ifdef PARANOID
1030  // Double check that the same number of coordinates (nsegments)
1031  // have been established for the boundary
1032  const unsigned nsegments_initial_coordinates =
1033  original_mesh_pt->boundary_segment_initial_coordinate(b).size();
1034 
1035  const unsigned nsegments_final_coordinates =
1036  original_mesh_pt->boundary_segment_final_coordinate(b).size();
1037 
1038  if (nsegments_initial_coordinates != nsegments_final_coordinates)
1039  {
1040  std::stringstream error_message;
1041  error_message
1042  << "The number of segments that present initial coordinates "
1043  << nsegments_initial_coordinates << " is different from "
1044  << "the\nnumber of segments that present final coordinates "
1045  << nsegments_final_coordinates << "\n\n";
1046  throw OomphLibError(error_message.str(),
1047  OOMPH_CURRENT_FUNCTION,
1048  OOMPH_EXCEPTION_LOCATION);
1049  } // if (nsegments_initial_coordinates!=nsegments_final_coordinates)
1050 
1051  // Also check that the number of segments found in the previous
1052  // mesh is the same as the number of segments found in this mesh
1053  if (nsegments_initial_coordinates != nsegments)
1054  {
1055  std::stringstream error_message;
1056  error_message << "Working with boundary (" << b
1057  << ").\n The number of initial and "
1058  << "final coordinates (" << nsegments_initial_coordinates
1059  << ") is different from\n"
1060  << "the number of found segments (" << nsegments
1061  << ").\n\n";
1062  throw OomphLibError(error_message.str(),
1063  OOMPH_CURRENT_FUNCTION,
1064  OOMPH_EXCEPTION_LOCATION);
1065  } // if (nsegments_initial_coordinates != nsegments)
1066 #endif
1067 
1068  // Create a backup for the data from the original mesh
1069  // Backup for the coordinates
1070  Vector<Vector<double>> original_mesh_segment_initial_coordinate(
1071  nsegments);
1072  Vector<Vector<double>> original_mesh_segment_final_coordinate(nsegments);
1073  // Backup for the zeta values
1074  Vector<double> original_mesh_segment_initial_zeta(nsegments);
1075  Vector<double> original_mesh_segment_final_zeta(nsegments);
1076  // Backup for the arclengths
1077  Vector<double> original_mesh_segment_initial_arclength(nsegments);
1078  Vector<double> original_mesh_segment_final_arclength(nsegments);
1079  // Do the backup
1080  for (unsigned is = 0; is < nsegments; is++)
1081  {
1082  original_mesh_segment_initial_coordinate[is].resize(2);
1083  original_mesh_segment_final_coordinate[is].resize(2);
1084  for (unsigned k = 0; k < 2; k++)
1085  {
1086  original_mesh_segment_initial_coordinate[is][k] =
1087  original_mesh_pt->boundary_segment_initial_coordinate(b)[is][k];
1088  original_mesh_segment_final_coordinate[is][k] =
1089  original_mesh_pt->boundary_segment_final_coordinate(b)[is][k];
1090  }
1091  // Check if the boudary has an associated GeomObject
1092  if (this->boundary_geom_object_pt(b) != 0)
1093  {
1094  original_mesh_segment_initial_zeta[is] =
1095  original_mesh_pt->boundary_segment_initial_zeta(b)[is];
1096  original_mesh_segment_final_zeta[is] =
1097  original_mesh_pt->boundary_segment_final_zeta(b)[is];
1098  }
1099  else
1100  {
1101  original_mesh_segment_initial_arclength[is] =
1102  original_mesh_pt->boundary_segment_initial_arclength(b)[is];
1103  original_mesh_segment_final_arclength[is] =
1104  original_mesh_pt->boundary_segment_final_arclength(b)[is];
1105  }
1106  } // for (is < nsegments)
1107 
1108  // Clear all the storage
1109  Boundary_segment_inverted[b].clear();
1110  Boundary_segment_initial_coordinate[b].clear();
1111  Boundary_segment_final_coordinate[b].clear();
1112 
1113  Boundary_segment_initial_zeta[b].clear();
1114  Boundary_segment_final_zeta[b].clear();
1115 
1116  Boundary_segment_initial_arclength[b].clear();
1117  Boundary_segment_final_arclength[b].clear();
1118 
1119  // Identify each segment in the processor with the ones created
1120  // by the original mesh
1121  // -----------------------------------------------------------------
1122  // Keep track of the already identified segments
1123  std::map<unsigned, bool> segment_done;
1124  for (unsigned is = 0; is < nsegments; is++)
1125  {
1126 #ifdef PARANOID
1127  // Flag to know if the segment was identified
1128  bool found_original_segment = false;
1129 #endif
1130 
1131  // Get the initial and final coordinates of the current segment
1132  Vector<double> current_seg_initial_coord(2);
1133  Vector<double> current_seg_final_coord(2);
1134 
1135  // Get access to the initial element on the segment
1136  FiniteElement* current_seg_initial_ele_pt =
1137  segment_sorted_ele_pt[is].front();
1138 
1139  // Number of nodes
1140  const unsigned nnod = current_seg_initial_ele_pt->nnode();
1141 
1142  // Get the first node of the current segment
1143  Node* current_seg_first_node_pt =
1144  current_seg_initial_ele_pt->node_pt(0);
1145  if (is_inverted[current_seg_initial_ele_pt])
1146  {
1147  current_seg_first_node_pt =
1148  current_seg_initial_ele_pt->node_pt(nnod - 1);
1149  }
1150 
1151  // Get access to the last element on the segment
1152  FiniteElement* current_seg_last_ele_pt =
1153  segment_sorted_ele_pt[is].back();
1154 
1155  // Get the last node of the current segment
1156  Node* current_seg_last_node_pt =
1157  current_seg_last_ele_pt->node_pt(nnod - 1);
1158  if (is_inverted[current_seg_last_ele_pt])
1159  {
1160  current_seg_last_node_pt = current_seg_last_ele_pt->node_pt(0);
1161  }
1162 
1163  // Get the coordinates for the first and last seg node
1164  for (unsigned i = 0; i < 2; i++)
1165  {
1166  current_seg_initial_coord[i] = current_seg_first_node_pt->x(i);
1167  current_seg_final_coord[i] = current_seg_last_node_pt->x(i);
1168  }
1169 
1170  // We have got the initial and final coordinates of the current
1171  // segment, compare those with the initial and final coordinates
1172  // of the original mesh segments to identify which segments is
1173  // which
1174  for (unsigned orig_s = 0; orig_s < nsegments; orig_s++)
1175  {
1176  if (!segment_done[orig_s])
1177  {
1178  // Get the coordinates to compare
1179  Vector<double> initial_coordinate =
1180  original_mesh_segment_initial_coordinate[orig_s];
1181  Vector<double> final_coordinate =
1182  original_mesh_segment_final_coordinate[orig_s];
1183 
1184  // Compute the distance initial(current)-initial(original)
1185  // coordinates
1186  double dist =
1187  ((current_seg_initial_coord[0] - initial_coordinate[0]) *
1188  (current_seg_initial_coord[0] - initial_coordinate[0])) +
1189  ((current_seg_initial_coord[1] - initial_coordinate[1]) *
1190  (current_seg_initial_coord[1] - initial_coordinate[1]));
1191  dist = sqrt(dist);
1192 
1193  // If the initial node is the same, check for the last node
1195  {
1196  // Compute the distance final(current)-final(original)
1197  // coordinates
1198  dist = ((current_seg_final_coord[0] - final_coordinate[0]) *
1199  (current_seg_final_coord[0] - final_coordinate[0])) +
1200  ((current_seg_final_coord[1] - final_coordinate[1]) *
1201  (current_seg_final_coord[1] - final_coordinate[1]));
1202  dist = sqrt(dist);
1203 
1204  // The final node is the same, we have identified the
1205  // segments
1207  {
1208  // Store the index that relates the previous index with the
1209  // current one
1210  current_segment_to_original_segment_index[is] = orig_s;
1211 
1212  // In this case the segment is not inverted
1213  Boundary_segment_inverted[b].push_back(0);
1214 
1215  // Copy the initial and final coordinates for each segment
1216  Boundary_segment_initial_coordinate[b].push_back(
1217  initial_coordinate);
1218  Boundary_segment_final_coordinate[b].push_back(
1219  final_coordinate);
1220 
1221  // Check if the boundary has an associated GeomObject
1222  if (this->boundary_geom_object_pt(b) != 0)
1223  {
1224  // Copy the initial zeta value for the segment
1225  Boundary_segment_initial_zeta[b].push_back(
1226  original_mesh_segment_initial_zeta[orig_s]);
1227  Boundary_segment_final_zeta[b].push_back(
1228  original_mesh_segment_final_zeta[orig_s]);
1229  }
1230  else
1231  {
1232  // Copy the initial and final arclength for each
1233  // segment
1234  Boundary_segment_initial_arclength[b].push_back(
1235  original_mesh_segment_initial_arclength[orig_s]);
1236  Boundary_segment_final_arclength[b].push_back(
1237  original_mesh_segment_final_arclength[orig_s]);
1238  }
1239  // Mark the segment as done
1240  segment_done[orig_s] = true;
1241 #ifdef PARANOID
1242  found_original_segment = true;
1243 #endif
1244  break;
1245  } // The final(current) node matched with the
1246  // final(original) node
1247  } // The initial(current) node matched with the
1248  // initial(original) node
1249  else
1250  {
1251  // Check the inverted case Compute the distance
1252  // initial(current)-final(original) coordinates
1253  double dist_inv =
1254  ((current_seg_initial_coord[0] - final_coordinate[0]) *
1255  (current_seg_initial_coord[0] - final_coordinate[0])) +
1256  ((current_seg_initial_coord[1] - final_coordinate[1]) *
1257  (current_seg_initial_coord[1] - final_coordinate[1]));
1258  dist_inv = sqrt(dist_inv);
1259 
1260  // If the initial node is the same as the final node of
1261  // the segment, check for the last node
1262  if (dist_inv <
1264  {
1265  // Compute the distance final(current)-initial(original)
1266  // coordinates
1267  dist_inv =
1268  ((current_seg_final_coord[0] - initial_coordinate[0]) *
1269  (current_seg_final_coord[0] - initial_coordinate[0])) +
1270  ((current_seg_final_coord[1] - initial_coordinate[1]) *
1271  (current_seg_final_coord[1] - initial_coordinate[1]));
1272  dist_inv = sqrt(dist_inv);
1273 
1274  // The final node is the same as the initial node, we
1275  // have identified the segments
1276  if (dist_inv <
1278  {
1279  // Store the index that related the previous index with the
1280  // current one
1281  current_segment_to_original_segment_index[is] = orig_s;
1282 
1283  // In this case the segment is inverted
1284  Boundary_segment_inverted[b].push_back(1);
1285 
1286  // Copy the initial and final coordinates for each segment
1287  Boundary_segment_initial_coordinate[b].push_back(
1288  initial_coordinate);
1289  Boundary_segment_final_coordinate[b].push_back(
1290  final_coordinate);
1291 
1292  // Check that the boudary has an associated GeomObject
1293  if (this->boundary_geom_object_pt(b) != 0)
1294  {
1295  // Copy the initial zeta value for the segments
1296  Boundary_segment_initial_zeta[b].push_back(
1297  original_mesh_segment_initial_zeta[orig_s]);
1298  Boundary_segment_final_zeta[b].push_back(
1299  original_mesh_segment_final_zeta[orig_s]);
1300  }
1301  else
1302  {
1303  // Copy the initial and final arclength for each segment
1304  Boundary_segment_initial_arclength[b].push_back(
1305  original_mesh_segment_initial_arclength[orig_s]);
1306  Boundary_segment_final_arclength[b].push_back(
1307  original_mesh_segment_final_arclength[orig_s]);
1308  }
1309  // Mark the segment as done
1310  segment_done[orig_s] = true;
1311 #ifdef PARANOID
1312  found_original_segment = true;
1313 #endif
1314  break;
1315  } // The final(current) node matched with the
1316  // initial(original) node
1317  } // The initial(current) node matched with the
1318  // final(original) node
1319  } // else (the first(current) node did not matched with the
1320  // first(original) node. Else do the inverted case
1321 
1322  } // (!segment_done[orig_s])
1323 
1324  } // (orig_s < nsegments)
1325 
1326 #ifdef PARANOID
1327  if (!found_original_segment)
1328  {
1329  std::stringstream error_message;
1330  error_message
1331  << "The (" << is << ")-th segment on the current segment was not\n"
1332  << "found when trying to identify it with the original mesh's\n"
1333  << "segment coordinates\n";
1334  throw OomphLibError(error_message.str(),
1335  OOMPH_CURRENT_FUNCTION,
1336  OOMPH_EXCEPTION_LOCATION);
1337  } // if (!found_original_segment)
1338 #endif
1339  } // for (is < nsegments)
1340 
1341  } // if (nsegments > 0)
1342 
1343  // -------------------------------------------------------------------
1344  // Fourth: The original mesh is different from the current mesh
1345  // (this). For boundaries with no geom object associated check if it
1346  // is required to reverse the zeta values. In order to reverse the
1347  // zeta values it is required to previously compute the arclength of
1348  // the segments and store the nodes in a container (set). NOTE that
1349  // the setup_boundary_coordinate() method is not called for
1350  // boundaries with NO GeomObject associated, so this is the LAST
1351  // CHANCE to do it
1352  // -------------------------------------------------------------------
1353  // The original mesh is the same as the current mesh (this). The
1354  // setup_boundary_method() will be called only for the boundaries
1355  // with NO GeomObject associated
1356  // -------------------------------------------------------------------
1357  if (this != original_mesh_pt)
1358  {
1359  // Get the boundary arclength
1360 
1361  // Get the initial and final zeta values for the boundary
1362  // (arclength) from the original mesh
1363  Vector<double> first_node_zeta_coordinate =
1364  original_mesh_pt->boundary_initial_zeta_coordinate(b);
1365  Vector<double> last_node_zeta_coordinate =
1366  original_mesh_pt->boundary_final_zeta_coordinate(b);
1367 
1368  // The boundary arclength is the maximum of the initial and final
1369  // zeta coordinate
1370  const double boundary_arclength =
1371  std::max(first_node_zeta_coordinate[0], last_node_zeta_coordinate[0]);
1372 
1373  for (unsigned is = 0; is < nsegments; is++)
1374  {
1375  // Here check if need to invert the elements and the boundary
1376  // coordinates for the segments in a boundary with no GeomObject
1377  // associated
1378  if (boundary_geom_object_pt(b) == 0)
1379  {
1380  // This case only applies for the initial and iterative mesh in
1381  // the adaptation process because the method
1382  // setup_boundary_coordinates() is called by the original mesh
1383  // for boundaries with no GeomObject associated
1384 
1385  // We are goind to check if it is necessary to invert the order
1386  // of the zeta values
1387 
1388  // Get the first and last node of the current segment and their
1389  // zeta values (arclength)
1390 
1391  // There is no need to check for nonhalo elements since the
1392  // container has only nonhalo face elements
1393 
1394  // Get access to the first element on the segment
1395  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
1396 
1397  // Number of nodes
1398  const unsigned nnod = first_ele_pt->nnode();
1399 
1400  // Get the first node of the current segment
1401  Node* first_node_pt = first_ele_pt->node_pt(0);
1402  if (is_inverted[first_ele_pt])
1403  {
1404  first_node_pt = first_ele_pt->node_pt(nnod - 1);
1405  }
1406 
1407  // Get access to the last element on the segment
1408  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
1409 
1410  // Get the last node of the current segment
1411  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
1412  if (is_inverted[last_ele_pt])
1413  {
1414  last_node_pt = last_ele_pt->node_pt(0);
1415  }
1416 
1417  // Get the zeta coordinates for the first and last node
1418  Vector<double> current_segment_initial_arclen(1);
1419  Vector<double> current_segment_final_arclen(1);
1420  // Is the segment in the current mesh (this) inverted?
1421  if (!Boundary_segment_inverted[b][is]) // Not inverted
1422  {
1423  first_node_pt->get_coordinates_on_boundary(
1424  b, current_segment_initial_arclen);
1425  last_node_pt->get_coordinates_on_boundary(
1426  b, current_segment_final_arclen);
1427  }
1428  else // Inverted
1429  {
1430  first_node_pt->get_coordinates_on_boundary(
1431  b, current_segment_final_arclen);
1432  last_node_pt->get_coordinates_on_boundary(
1433  b, current_segment_initial_arclen);
1434  }
1435 
1436  // Once the zeta values have been obtained check if they are set
1437  // in increasing or decreasing order
1438 
1439  // Flag to state that the values in the segment are in increasing
1440  // order
1441  bool increasing_order = false;
1442 
1443  // If the initial zeta value is smaller than the final zeta
1444  // value then they are in increasing order
1445  if (current_segment_initial_arclen[0] <
1446  current_segment_final_arclen[0])
1447  {
1448  increasing_order = true;
1449  }
1450  // If the initial zeta value is greater than the initial zeta
1451  // value then they are in decreasing order
1452  else if (current_segment_initial_arclen[0] >
1453  current_segment_final_arclen[0])
1454  {
1455  increasing_order = false;
1456  }
1457 #ifdef PARANOID
1458  else
1459  {
1460  std::stringstream error_message;
1461  error_message
1462  << "It was not possible to identify if the zeta values on "
1463  << "boundary (" << b << ")\nand segment (" << is
1464  << ") should go in "
1465  << "increasing or decreasing order.\n--- New mesh ---\n"
1466  << "Current segment initial arclength: ("
1467  << current_segment_initial_arclen[0] << ")\n"
1468  << "First node coordinates: (" << first_node_pt->x(0) << ", "
1469  << first_node_pt->x(1) << ")\n"
1470  << "Current segment final arclength: ("
1471  << current_segment_final_arclen[0] << ")\n"
1472  << "Last node coordinates: (" << last_node_pt->x(0) << ", "
1473  << last_node_pt->x(1) << ")\n"
1474  << "Current segment arclength: (" << segment_arclength[is]
1475  << ")\n";
1476  throw OomphLibError(error_message.str(),
1477  OOMPH_CURRENT_FUNCTION,
1478  OOMPH_EXCEPTION_LOCATION);
1479  }
1480 #endif
1481 
1482  // Now get the original initial and final arclengths and check
1483  // if they are in increasing or decreasing order
1484  const unsigned prev_s = current_segment_to_original_segment_index[is];
1485  const double original_segment_initial_arclength =
1486  original_mesh_pt->boundary_segment_initial_arclength(b)[prev_s];
1487  const double original_segment_final_arclength =
1488  original_mesh_pt->boundary_segment_final_arclength(b)[prev_s];
1489 
1490  // Flag to check if the values go in increasing or decreasing
1491  // order in the original mesh segment
1492  bool original_increasing_order = false;
1493 
1494  // Now check if the arclengths on the original mesh go in
1495  // increase or decrease order, this is also used to choose the
1496  // starting value to map the values in the current segment
1497  double starting_arclength = 0.0;
1498  if (original_segment_final_arclength >
1499  original_segment_initial_arclength)
1500  {
1501  // ... in increasing order in the original mesh ...
1502  original_increasing_order = true;
1503  // Select the starting arclength
1504  starting_arclength = original_segment_initial_arclength;
1505  }
1506  else if (original_segment_final_arclength <
1507  original_segment_initial_arclength)
1508  {
1509  // ... in decreasing order in the original mesh ...
1510  original_increasing_order = false;
1511  // Select the starting arclength
1512  starting_arclength = original_segment_final_arclength;
1513  }
1514 #ifdef PARANOID
1515  else
1516  {
1517  std::stringstream error_message;
1518  error_message
1519  << "It was not possible to identify if the zeta values on "
1520  << "boundary (" << b << ")\nand segment (" << is
1521  << ") should go in "
1522  << "increasing or decreasing order.\n--- Original mesh ---\n"
1523  << "Original segment initial arclength: ("
1524  << original_segment_initial_arclength << ")\n"
1525  << "Original segment final arclength: ("
1526  << original_segment_final_arclength << ")\n";
1527  throw OomphLibError(error_message.str(),
1528  OOMPH_CURRENT_FUNCTION,
1529  OOMPH_EXCEPTION_LOCATION);
1530  }
1531 #endif
1532 
1533  // Now scale the zeta values based considering if the zeta
1534  // values from the current mesh (this) go in the same order as
1535  // in the original mesh
1536  if (increasing_order && original_increasing_order)
1537  {
1538  // Current seg
1539  // |------|
1540  // 0 ---- 1
1541  //
1542  // Is mapped to the new values
1543  // |------|
1544  // a ---- b
1545  // a = original_segment_initial_arclength
1546  // b = original_segment_final_arclength
1547  // s = starting_arclength
1548  // The mapping is given by
1549  // new_z = s + z_old * (b - a)
1550 
1551  // Get the nodes associated to the segment
1552  std::set<Node*> seg_nodes_pt = segment_all_nodes_pt[is];
1553  // Go through all the nodes in the segment an change their
1554  // zeta values
1555  for (std::set<Node*>::iterator it = seg_nodes_pt.begin();
1556  it != seg_nodes_pt.end();
1557  it++)
1558  {
1559  // Storing for the zeta value
1560  Vector<double> zeta(1);
1561  // Get each node
1562  Node* nod_pt = (*it);
1563  // Get the zeta value of the current node
1564  nod_pt->get_coordinates_on_boundary(b, zeta);
1565  // ... and re-assign it
1566  const double temp =
1567  starting_arclength + (zeta[0] * segment_arclength[is]);
1568  // The zeta value
1569  zeta[0] = temp / boundary_arclength;
1570  // Correct
1571  if (std::fabs(zeta[0] - 1.0) < 1.0e-14)
1572  {
1573  zeta[0] = 1.0;
1574  }
1575  else if (std::fabs(zeta[0]) < 1.0e-14)
1576  {
1577  zeta[0] = 0.0;
1578  }
1579 
1580  // Set the new value
1581  nod_pt->set_coordinates_on_boundary(b, zeta);
1582  } // Go through all the nodes
1583  } // if (increasing_order && original_increasing_order)
1584  else if (!increasing_order && original_increasing_order)
1585  {
1586  // Current seg
1587  // |------|
1588  // 1 ---- 0
1589  //
1590  // Is mapped to the new values
1591  // |------|
1592  // a ---- b
1593  // a = original_segment_initial_arclength
1594  // b = original_segment_final_arclength
1595  // s = starting_arclength
1596  // The mapping is given by
1597  // new_z = s + (1.0 - z_old) * (b - a)
1598 
1599  // Get the nodes associated to the segment
1600  std::set<Node*> seg_nodes_pt = segment_all_nodes_pt[is];
1601  // Go through all the nodes in the segment an change their
1602  // zeta values
1603  for (std::set<Node*>::iterator it = seg_nodes_pt.begin();
1604  it != seg_nodes_pt.end();
1605  it++)
1606  {
1607  // Storing for the zeta value
1608  Vector<double> zeta(1);
1609  // Get each node
1610  Node* nod_pt = (*it);
1611  // Get the zeta value of the current node
1612  nod_pt->get_coordinates_on_boundary(b, zeta);
1613  // ... and re-assign it
1614  const double temp =
1615  starting_arclength + ((1.0 - zeta[0]) * segment_arclength[is]);
1616  // The zeta value
1617  zeta[0] = temp / boundary_arclength;
1618  // Correct
1619  if (std::fabs(zeta[0] - 1.0) < 1.0e-14)
1620  {
1621  zeta[0] = 1.0;
1622  }
1623  else if (std::fabs(zeta[0]) < 1.0e-14)
1624  {
1625  zeta[0] = 0.0;
1626  }
1627  // Set the new value
1628  nod_pt->set_coordinates_on_boundary(b, zeta);
1629  } // Go through all the nodes
1630  } // else if (!increasing_order && original_increasing_order)
1631  else if (increasing_order && !original_increasing_order)
1632  {
1633  // Current seg
1634  // |------|
1635  // 0 ---- 1
1636  //
1637  // Is mapped to the new values
1638  // |------|
1639  // b ---- a
1640  // a = original_segment_initial_arclength
1641  // b = original_segment_final_arclength
1642  // s = starting_arclength
1643  // The mapping is given by
1644  // new_z = s + (1.0 - z_old) * |(b - a)|
1645 
1646  // Get the nodes associated to the segment
1647  std::set<Node*> seg_nodes_pt = segment_all_nodes_pt[is];
1648  // Go through all the nodes in the segment an change their
1649  // zeta values
1650  for (std::set<Node*>::iterator it = seg_nodes_pt.begin();
1651  it != seg_nodes_pt.end();
1652  it++)
1653  {
1654  // Storing for the zeta value
1655  Vector<double> zeta(1);
1656  // Get each node
1657  Node* nod_pt = (*it);
1658  // Get the zeta value of the current node
1659  nod_pt->get_coordinates_on_boundary(b, zeta);
1660  // ... and re-assign it
1661  const double temp =
1662  starting_arclength + ((1.0 - zeta[0]) * segment_arclength[is]);
1663  // The zeta value
1664  zeta[0] = temp / boundary_arclength;
1665  // Correct
1666  if (std::fabs(zeta[0] - 1.0) < 1.0e-14)
1667  {
1668  zeta[0] = 1.0;
1669  }
1670  else if (std::fabs(zeta[0]) < 1.0e-14)
1671  {
1672  zeta[0] = 0.0;
1673  }
1674  // Set the new value
1675  nod_pt->set_coordinates_on_boundary(b, zeta);
1676  } // Go through all the nodes
1677  } // else if (increasing_order && !original_increasing_order)
1678  else if (!increasing_order && !original_increasing_order)
1679  {
1680  // Current seg
1681  // |------|
1682  // 0 ---- 1
1683  //
1684  // Is mapped to the new values
1685  // |------|
1686  // a ---- b
1687  // a = original_segment_initial_arclength
1688  // b = original_segment_final_arclength
1689  // s = starting_arclength
1690  // The mapping is given by
1691  // new_z = s + z_old * |(b - a)|
1692 
1693  // Get the nodes associated to the segment
1694  std::set<Node*> seg_nodes_pt = segment_all_nodes_pt[is];
1695  // Go through all the nodes in the segment an change their
1696  // zeta values
1697  for (std::set<Node*>::iterator it = seg_nodes_pt.begin();
1698  it != seg_nodes_pt.end();
1699  it++)
1700  {
1701  // Storing for the zeta value
1702  Vector<double> zeta(1);
1703  // Get each node
1704  Node* nod_pt = (*it);
1705  // Get the zeta value of the current node
1706  nod_pt->get_coordinates_on_boundary(b, zeta);
1707  // ... and re-assign it
1708  const double temp =
1709  starting_arclength + (zeta[0] * segment_arclength[is]);
1710  // The zeta value
1711  zeta[0] = temp / boundary_arclength;
1712  // Correct
1713  if (std::fabs(zeta[0] - 1.0) < 1.0e-14)
1714  {
1715  zeta[0] = 1.0;
1716  }
1717  else if (std::fabs(zeta[0]) < 1.0e-14)
1718  {
1719  zeta[0] = 0.0;
1720  }
1721  // Set the new value
1722  nod_pt->set_coordinates_on_boundary(b, zeta);
1723  } // Go through all the nodes
1724  } // else if (!increasing_order && !original_increasing_order)
1725 
1726 #ifdef PARANOID
1727  // Verify that the z values of the first and last node are not
1728  // out of the range [0,1]
1729  for (std::list<FiniteElement*>::iterator it_list =
1730  segment_sorted_ele_pt[is].begin();
1731  it_list != segment_sorted_ele_pt[is].end();
1732  it_list++)
1733  {
1734  // Number of nodes in the segment
1735  const unsigned nnod = (*it_list)->nnode();
1736 
1737  // Get the first node of the current segment
1738  Node* first_node_pt = (*it_list)->node_pt(0);
1739  if (is_inverted[(*it_list)])
1740  {
1741  first_node_pt = (*it_list)->node_pt(nnod - 1);
1742  }
1743 
1744  // Get the last node of the current segment
1745  Node* last_node_pt = (*it_list)->node_pt(nnod - 1);
1746  if (is_inverted[(*it_list)])
1747  {
1748  last_node_pt = (*it_list)->node_pt(0);
1749  }
1750 
1751  // The z value for the first node
1752  Vector<double> zeta(1);
1753  first_node_pt->get_coordinates_on_boundary(b, zeta);
1754  if (zeta[0] < 0.0 || zeta[0] > 1.0)
1755  {
1756  std::ostringstream error_message;
1757  error_message
1758  << "The boundary coordinate of the first node on boundary ("
1759  << b << ")\nand segment (" << is << ") is out of the "
1760  << "allowed values [0,1]\n"
1761  << "The node boundary coordinate: (" << zeta[0] << ")\n"
1762  << "The vertex coordinates are: (" << first_node_pt->x(0)
1763  << ", " << first_node_pt->x(1) << ")\n";
1764  throw OomphLibError(error_message.str(),
1765  OOMPH_CURRENT_FUNCTION,
1766  OOMPH_EXCEPTION_LOCATION);
1767  }
1768 
1769  // The z value for the last node
1770  last_node_pt->get_coordinates_on_boundary(b, zeta);
1771  if (zeta[0] < 0.0 || zeta[0] > 1.0)
1772  {
1773  std::ostringstream error_message;
1774  error_message
1775  << "The boundary coordinate of the last node on boundary (" << b
1776  << ")\nand segment (" << is << ") is out of the "
1777  << "allowed values [0,1]\n"
1778  << "The node boundary coordinate: (" << zeta[0] << ")\n"
1779  << "The vertex coordinates are: (" << last_node_pt->x(0) << ", "
1780  << last_node_pt->x(1) << ")\n";
1781  throw OomphLibError(error_message.str(),
1782  OOMPH_CURRENT_FUNCTION,
1783  OOMPH_EXCEPTION_LOCATION);
1784  }
1785  }
1786 #endif // #ifdef PARANOID
1787 
1788  } // if (boundary_geom_object_pt(b)==0)
1789 
1790  } // for (is < nsegments)
1791 
1792  } // if (this != original_mesh_pt)
1793 
1794  // ------------------------------------------------------------------
1795  // Copy the corrected (possible reversed) info. to the containers of
1796  // the current mesh
1797  // ------------------------------------------------------------------
1798  // Check if there are segments of b boundary in this processor
1799  if (nsegments > 0)
1800  {
1801  // Copy the initial and final coordinates
1802  Boundary_initial_coordinate[b] =
1803  original_mesh_pt->boundary_initial_coordinate(b);
1804 
1805  Boundary_final_coordinate[b] =
1806  original_mesh_pt->boundary_final_coordinate(b);
1807 
1808  // The initial and final zeta coordinates (In case of a geometric
1809  // object those are the limits of the geom object)
1810  Boundary_initial_zeta_coordinate[b] =
1811  original_mesh_pt->boundary_initial_zeta_coordinate(b);
1812 
1813  Boundary_final_zeta_coordinate[b] =
1814  original_mesh_pt->boundary_final_zeta_coordinate(b);
1815 
1816  } // if (nsegments > 0)
1817 
1818  // Set the flag to indicate that the zeta values have been assigned
1819  // for the current boundary
1820  Assigned_segments_initial_zeta_values[b] = true;
1821 
1822  // Clean all the created face elements
1823  for (unsigned i = 0; i < nele; i++)
1824  {
1825  delete face_el_pt[i];
1826  face_el_pt[i] = 0;
1827  }
1828  }
1829 
1830  //======================================================================
1831  /// Compute the boundary segments connectivity for those
1832  /// boundaries that were splited during the distribution process
1833  /// and also the initial zeta values for each segment (the initial
1834  /// and final boundary nodes coordinates)
1835  //======================================================================
1836  template<class ELEMENT>
1839  const unsigned& b)
1840  {
1841  // ------------------------------------------------------------------
1842  // First: Get the face elements associated with the current boundary
1843  // ------------------------------------------------------------------
1844 
1845  // Get the communicator of the mesh
1846  OomphCommunicator* comm_pt = this->communicator_pt();
1847 
1848  // Get the number of processors
1849  const unsigned nproc = comm_pt->nproc();
1850  // Get the rank of the current processor
1851  const unsigned my_rank = comm_pt->my_rank();
1852 
1853  // Temporary storage for face elements
1854  Vector<FiniteElement*> all_face_ele_pt;
1855 
1856  // Flag to know whether we are working with an internal open curve
1857  // and then re-assign the initial and final zeta coordinates for
1858  // each segment (only used when the mesh is distributed)
1859  bool is_internal_boundary = false;
1860 
1861  // map to associate the face element to the bulk element, necessary
1862  // to attach halo face elements at both sides of each found segment
1863  std::map<FiniteElement*, FiniteElement*> face_to_bulk_element_pt;
1864 
1865  // Select the boundary face elements, using the criteria of highest
1866  // processor in charge and bottom-left element
1867  select_boundary_face_elements(
1868  all_face_ele_pt, b, is_internal_boundary, face_to_bulk_element_pt);
1869 
1870  // Get the number of face elements
1871  const unsigned n_all_face_ele = all_face_ele_pt.size();
1872 
1873  // ----------------------------------------------------------------
1874  // Second: Sort the face elements, only consider nonhalo elements
1875  // ----------------------------------------------------------------
1876 
1877  // A flag vector to mark those face elements that are considered as
1878  // halo in the current processor
1879  std::vector<bool> is_halo_face_element(n_all_face_ele, false);
1880 
1881  // Count the total number of non halo face elements
1882  unsigned nnon_halo_face_elements = 0;
1883 
1884  // Only mark the face elements as halo if the mesh is marked as
1885  // distributed
1886  for (unsigned ie = 0; ie < n_all_face_ele; ie++)
1887  {
1888  FiniteElement* face_ele_pt = all_face_ele_pt[ie];
1889  // Get the bulk element
1890  FiniteElement* tmp_bulk_ele_pt = face_to_bulk_element_pt[face_ele_pt];
1891  // Check if the bulk element is halo
1892  if (!tmp_bulk_ele_pt->is_halo())
1893  {
1894  // Set the flag for non halo element
1895  is_halo_face_element[ie] = false;
1896  // Increase the non halo elements counter
1897  nnon_halo_face_elements++;
1898  }
1899  else
1900  {
1901  // Mark the face element as halo
1902  is_halo_face_element[ie] = true;
1903  }
1904 
1905  } // for (ie < n_ele)
1906 
1907  // Get the total number of halo face elements
1908  const unsigned nhalo_face_element =
1909  n_all_face_ele - nnon_halo_face_elements;
1910 
1911  // The vector of list to store the "segments" that compound the
1912  // boundary (segments may appear only in a distributed mesh)
1913  Vector<std::list<FiniteElement*>> segment_sorted_ele_pt;
1914 
1915  // Number of already sorted face elements (only nonhalo elements for
1916  // a distributed mesh)
1917  unsigned nsorted_face_elements = 0;
1918 
1919  // Keep track of who's done (this apply to nonhalo only, remember we
1920  // are only working with halo elements)
1921  std::map<FiniteElement*, bool> done_el;
1922 
1923  // Keep track of which element is inverted (in distributed mesh the
1924  // elements may be inverted with respect to the segment they belong)
1925  std::map<FiniteElement*, bool> is_inverted;
1926 
1927  // Iterate until all possible segments have been created
1928  while (nsorted_face_elements < nnon_halo_face_elements)
1929  {
1930  // The ordered list of face elements (in a distributed mesh a
1931  // collection of contiguous face elements define a segment)
1932  std::list<FiniteElement*> sorted_el_pt;
1933  sorted_el_pt.clear();
1934 
1935 #ifdef PARANOID
1936  // Select an initial element for the segment (the first not done
1937  // nonhalo element)
1938  bool found_initial_face_element = false;
1939 #endif
1940 
1941  FiniteElement* ele_face_pt = 0;
1942 
1943  unsigned iface = 0;
1944  for (iface = 0; iface < n_all_face_ele; iface++)
1945  {
1946  if (!is_halo_face_element[iface])
1947  {
1948  ele_face_pt = all_face_ele_pt[iface];
1949  // If not done then take it as initial face element
1950  if (!done_el[ele_face_pt])
1951  {
1952 #ifdef PARANOID
1953  found_initial_face_element = true;
1954 #endif
1955  nsorted_face_elements++;
1956  iface++; // The next element number
1957  sorted_el_pt.push_back(ele_face_pt);
1958  // Mark as done
1959  done_el[ele_face_pt] = true;
1960  break;
1961  }
1962  }
1963  } // for (iface < nele)
1964 
1965 #ifdef PARANOID
1966  if (!found_initial_face_element)
1967  {
1968  std::ostringstream error_message;
1969  error_message
1970  << "Could not find an initial face element for the current segment\n";
1971  // << "----- Possible memory leak -----\n";
1972  throw OomphLibError(error_message.str(),
1973  OOMPH_CURRENT_FUNCTION,
1974  OOMPH_EXCEPTION_LOCATION);
1975  }
1976 #endif
1977 
1978  // Number of nodes
1979  const unsigned nnod = ele_face_pt->nnode();
1980 
1981  // Left and rightmost nodes (the left and right nodes of the
1982  // current face element)
1983  Node* left_node_pt = ele_face_pt->node_pt(0);
1984  Node* right_node_pt = ele_face_pt->node_pt(nnod - 1);
1985 
1986  // Continue iterating if a new face element has been added to the
1987  // list
1988  bool face_element_added = false;
1989 
1990  // While a new face element has been added to the set of sorted
1991  // face elements then re-iterate
1992  do
1993  {
1994  // Start from the next face element since we have already added
1995  // the previous one as the initial face element (any previous
1996  // face element had to be added on previous iterations)
1997  for (unsigned iiface = iface; iiface < n_all_face_ele; iiface++)
1998  {
1999  // Re-start flag
2000  face_element_added = false;
2001 
2002  // Get the candidate element
2003  ele_face_pt = all_face_ele_pt[iiface];
2004 
2005  // Check that the candidate element has not been done and is
2006  // not a halo element
2007  if (!(done_el[ele_face_pt] || is_halo_face_element[iiface]))
2008  {
2009  // Get the left and right nodes of the current element
2010  Node* local_left_node_pt = ele_face_pt->node_pt(0);
2011  Node* local_right_node_pt = ele_face_pt->node_pt(nnod - 1);
2012 
2013  // New element fits at the left of segment and is not inverted
2014  if (left_node_pt == local_right_node_pt)
2015  {
2016  left_node_pt = local_left_node_pt;
2017  sorted_el_pt.push_front(ele_face_pt);
2018  is_inverted[ele_face_pt] = false;
2019  face_element_added = true;
2020  }
2021  // New element fits at the left of segment and is inverted
2022  else if (left_node_pt == local_left_node_pt)
2023  {
2024  left_node_pt = local_right_node_pt;
2025  sorted_el_pt.push_front(ele_face_pt);
2026  is_inverted[ele_face_pt] = true;
2027  face_element_added = true;
2028  }
2029  // New element fits on the right of segment and is not inverted
2030  else if (right_node_pt == local_left_node_pt)
2031  {
2032  right_node_pt = local_right_node_pt;
2033  sorted_el_pt.push_back(ele_face_pt);
2034  is_inverted[ele_face_pt] = false;
2035  face_element_added = true;
2036  }
2037  // New element fits on the right of segment and is inverted
2038  else if (right_node_pt == local_right_node_pt)
2039  {
2040  right_node_pt = local_left_node_pt;
2041  sorted_el_pt.push_back(ele_face_pt);
2042  is_inverted[ele_face_pt] = true;
2043  face_element_added = true;
2044  }
2045 
2046  if (face_element_added)
2047  {
2048  done_el[ele_face_pt] = true;
2049  nsorted_face_elements++;
2050  break;
2051  }
2052 
2053  } // if (!(done_el[ele_face_pt] || is_halo_face_element[iiface]))
2054  } // for (iiface<nnon_halo_face_element)
2055  } while (face_element_added &&
2056  (nsorted_face_elements < nnon_halo_face_elements));
2057 
2058  // Store the created segment in the vector of segments
2059  segment_sorted_ele_pt.push_back(sorted_el_pt);
2060 
2061  } // while(nsorted_face_elements < nnon_halo_face_elements);
2062 
2063  // -----------------------------------------------------------------
2064  // Third: We have the face elements sorted (in segments), now assign
2065  // boundary coordinates to the nodes in the segments, this is the
2066  // LOCAL boundary coordinate and further communication is needed to
2067  // compute the GLOBAL boundary coordinates
2068  // -----------------------------------------------------------------
2069 
2070  // Vector of sets that stores the nodes of each segment based on a
2071  // lexicographically order starting from the bottom left node of
2072  // each segment
2073  Vector<std::set<Node*>> segment_all_nodes_pt;
2074 
2075  // The number of segments in this processor
2076  const unsigned nsegments = segment_sorted_ele_pt.size();
2077  // DEBP(nsegments);
2078 
2079 #ifdef PARANOID
2080  if (nnon_halo_face_elements > 0 && nsegments == 0)
2081  {
2082  std::ostringstream error_message;
2083  error_message
2084  << "The number of segments is zero, but the number of nonhalo\n"
2085  << "elements is: (" << nnon_halo_face_elements << ")\n";
2086  throw OomphLibError(
2087  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2088  } // if (nnon_halo_face_elements > 0 && nsegments == 0)
2089 #endif
2090 
2091  // The arclength of each segment in the current processor
2092  Vector<double> segment_arclength(nsegments);
2093 
2094  // The number of vertices of each segment
2095  Vector<unsigned> nvertices_per_segment(nsegments);
2096 
2097  // The initial zeta for the segment
2098  Vector<double> initial_zeta_segment(nsegments);
2099 
2100  // The final zeta for the segment
2101  Vector<double> final_zeta_segment(nsegments);
2102 
2103  // Go through all the segments and compute its ARCLENGTH (if the
2104  // boundary has a GeomObject associated then assign the initial and
2105  // final zeta values for the segment)
2106  for (unsigned is = 0; is < nsegments; is++)
2107  {
2108 #ifdef PARANOID
2109  if (segment_sorted_ele_pt[is].size() == 0)
2110  {
2111  std::ostringstream error_message;
2112  error_message << "The (" << is << ")-th segment has no elements\n";
2113  throw OomphLibError(error_message.str(),
2114  OOMPH_CURRENT_FUNCTION,
2115  OOMPH_EXCEPTION_LOCATION);
2116  } // if (segment_sorted_ele_pt[is].size() == 0)
2117 #endif
2118 
2119  // Get access to the first element on the segment
2120  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
2121 
2122  // Number of nodes
2123  const unsigned nnod = first_ele_pt->nnode();
2124 
2125  // Get the first node of the current segment
2126  Node* first_node_pt = first_ele_pt->node_pt(0);
2127  if (is_inverted[first_ele_pt])
2128  {
2129  first_node_pt = first_ele_pt->node_pt(nnod - 1);
2130  }
2131 
2132  // Coordinates of left node
2133  double x_left = first_node_pt->x(0);
2134  double y_left = first_node_pt->x(1);
2135 
2136  // Initialise boundary coordinate (local boundary coordinate for
2137  // boundaries with more than one segment)
2138  Vector<double> zeta(1, 0.0);
2139 
2140  // If we have associated a GeomObject then it is not necessary to
2141  // compute the arclength, only read the values from the nodes at
2142  // the edges and set the initial and final zeta segment values
2143  if (this->boundary_geom_object_pt(b) != 0)
2144  {
2145  // Get the initial node coordinate
2146  first_node_pt->get_coordinates_on_boundary(b, zeta);
2147  // Set the initial zeta segment value
2148  initial_zeta_segment[is] = zeta[0];
2149 
2150  // Get access to the last element on the segment
2151  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
2152 
2153  // Get the last node of the current segment
2154  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
2155  if (is_inverted[last_ele_pt])
2156  {
2157  last_node_pt = last_ele_pt->node_pt(0);
2158  }
2159 
2160  // Get the final node coordinate
2161  last_node_pt->get_coordinates_on_boundary(b, zeta);
2162  // Set the final zeta segment value
2163  final_zeta_segment[is] = zeta[0];
2164  }
2165 
2166  // Sort the nodes in the segment (lexicographically bottom left
2167  // node)
2168  std::set<Node*> local_nodes_pt;
2169  // Insert the first node
2170  local_nodes_pt.insert(first_node_pt);
2171 
2172  // Now loop over nodes in order and increase the ARCLENGTH
2173  for (std::list<FiniteElement*>::iterator it =
2174  segment_sorted_ele_pt[is].begin();
2175  it != segment_sorted_ele_pt[is].end();
2176  it++)
2177  {
2178  // Get the pointer to the element
2179  FiniteElement* el_pt = (*it);
2180 
2181  // Start node and increment
2182  unsigned k_nod = 1;
2183  int nod_diff = 1;
2184  // Access nodes in reverse?
2185  if (is_inverted[el_pt])
2186  {
2187  k_nod = nnod - 2;
2188  nod_diff = -1;
2189  }
2190 
2191  // Loop over nodes in the face element
2192  for (unsigned j = 1; j < nnod; j++)
2193  {
2194  Node* nod_pt = el_pt->node_pt(k_nod);
2195  k_nod += nod_diff;
2196 
2197  // Coordinates of right node
2198  double x_right = nod_pt->x(0);
2199  double y_right = nod_pt->x(1);
2200 
2201  // Increment boundary coordinate (the arclength)
2202  zeta[0] += sqrt((x_right - x_left) * (x_right - x_left) +
2203  (y_right - y_left) * (y_right - y_left));
2204 
2205  // When we have a GeomObject associated to the boundary we already
2206  // know the zeta values for the nodes, there is no need to compute
2207  // the arclength
2208  // if (this->boundary_geom_object_pt(b)==0)
2209  // {
2210  // // Set boundary coordinate
2211  // // nod_pt->set_coordinates_on_boundary(b, zeta);
2212  // }
2213 
2214  // Increment reference coordinate
2215  x_left = x_right;
2216  y_left = y_right;
2217 
2218  // Get lexicographically bottom left node but only
2219  // use vertex nodes as candidates
2220  local_nodes_pt.insert(nod_pt);
2221 
2222  } // for (j < nnod)
2223 
2224  } // iterator over the elements in the segment
2225 
2226  // Info. to be passed to other processors
2227  // The initial arclength for the segment that goes after this depends
2228  // on the current segment arclength
2229  segment_arclength[is] = zeta[0];
2230 
2231  // Info. to be passed to the other processors
2232  // The initial vertex number for the segment that goes after this
2233  // depends on the current segment vertices number
2234  nvertices_per_segment[is] = local_nodes_pt.size();
2235 
2236  // Add the nodes for the corresponding segment in the container
2237  segment_all_nodes_pt.push_back(local_nodes_pt);
2238 
2239  // The attaching of the halo elements at both sides of the segments is
2240  // performed only if segments connectivity needs to be computed
2241 
2242  } // for (is < nsegments)
2243 
2244  // Container to store the number of vertices before each segment,
2245  // initialise to zero in case we have a non distributed boundary
2246  Vector<unsigned> nvertices_before_segment(nsegments, 0);
2247 
2248  // Store the initial arclength for each segment of boundary in the
2249  // current processor, initalise to zero in case we have a non
2250  // distributed boundary
2251  Vector<double> initial_segment_arclength(nsegments, 0.0);
2252 
2253  // Info. to be passed to other processors
2254  // If the boundary is distributed we need to know which processors does
2255  // have the initial and final segments, this helps to get the first and
2256  // last nodes coordinates (info. used to scale the bound coordinates)
2257 
2258  // Processors with the initial and final segment
2259  unsigned proc_with_initial_seg = 0;
2260  unsigned proc_with_final_seg = 0;
2261 
2262  // ... and the index of those segments (only of interest in the
2263  // processors that have the initial and final segments)
2264  unsigned initial_segment = 0;
2265  unsigned final_segment = 0;
2266 
2267  // Each segment needs to know whether it has to be inverted or not
2268  // Store whether a segment needs to be inverted or not
2269  Vector<unsigned> segment_inverted(nsegments);
2270 
2271  // Before attaching the halo elements create a copy of the data
2272  // structure without halo elements
2273  Vector<std::list<FiniteElement*>> segment_sorted_nonhalo_ele_pt(nsegments);
2274  for (unsigned is = 0; is < nsegments; is++)
2275  {
2276  for (std::list<FiniteElement*>::iterator it_seg =
2277  segment_sorted_ele_pt[is].begin();
2278  it_seg != segment_sorted_ele_pt[is].end();
2279  it_seg++)
2280  {
2281  segment_sorted_nonhalo_ele_pt[is].push_back((*it_seg));
2282  }
2283 
2284  } // for (is < nsegments)
2285 
2286  // --------------------------------------------------------------
2287  // Attach the halo elements at both sides of the segments
2288  for (unsigned is = 0; is < nsegments; is++)
2289  {
2290  // Get access to the first element on the segment
2291  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
2292 
2293  // Number of nodes
2294  const unsigned nnod = first_ele_pt->nnode();
2295 
2296  // Get the first node of the current segment
2297  Node* first_node_pt = first_ele_pt->node_pt(0);
2298  if (is_inverted[first_ele_pt])
2299  {
2300  first_node_pt = first_ele_pt->node_pt(nnod - 1);
2301  }
2302 
2303  // Get access to the last element on the segment
2304  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
2305 
2306  // Get the last node of the current segment
2307  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
2308  if (is_inverted[last_ele_pt])
2309  {
2310  last_node_pt = last_ele_pt->node_pt(0);
2311  }
2312 
2313  // -----------------------------------------------------------------
2314  // Fourth: Now attach the halo elements to the left and right side
2315  // of each segment
2316  // -----------------------------------------------------------------
2317  bool attached_left_halo = false;
2318  bool attached_right_halo = false;
2319  if (nhalo_face_element > 0)
2320  {
2321  for (unsigned iiface = 0; iiface < n_all_face_ele; iiface++)
2322  {
2323  // Get the candidate element
2324  FiniteElement* halo_face_ele_pt = all_face_ele_pt[iiface];
2325 
2326  // Check that the element is a halo face element, we do not check
2327  // if the element has been already done since the halo elements
2328  // may be connected to more than one segment (2 at most), to the
2329  // left and right of different segments
2330  //
2331  // Segment k Halo Segment r
2332  // |---|---|---| |xxx| |---|---|---|
2333  //
2334  // Segment k Halo Segment r
2335  // |---|---|---|xxx|---|---|---|
2336  //
2337  if (is_halo_face_element[iiface])
2338  {
2339  // Get its left and right nodes
2340  Node* left_node_pt = halo_face_ele_pt->node_pt(0);
2341  Node* right_node_pt = halo_face_ele_pt->node_pt(nnod - 1);
2342  // The halo element fits to the left of segment
2343  if (!attached_left_halo && (first_node_pt == right_node_pt ||
2344  first_node_pt == left_node_pt))
2345  {
2346  // Add the halo element to the left of the segment
2347  segment_sorted_ele_pt[is].push_front(halo_face_ele_pt);
2348 
2349  // Once a halo face element has been added to the left
2350  // mark as found halo to the left
2351  attached_left_halo = true;
2352  }
2353  // The halo element fits to the right of the segment
2354  else if (!attached_right_halo && (last_node_pt == left_node_pt ||
2355  last_node_pt == right_node_pt))
2356  {
2357  // Add the halo element to the right of the segment
2358  segment_sorted_ele_pt[is].push_back(halo_face_ele_pt);
2359  // Once a halo face element has been added to the right
2360  // mark as found halo to the right
2361  attached_right_halo = true;
2362  }
2363  // If we have already found elements to left and right then
2364  // break the loop
2365  if (attached_left_halo && attached_right_halo)
2366  {
2367  break;
2368  }
2369 
2370  } // if (is_halo_face_element[iiface])
2371 
2372  } // for (iiface < nel)
2373 
2374  } // if (nhalo_face_element > 0)
2375 
2376  } // for (is < nsegments)
2377 
2378  // The segments now have local coordinates assigned and halo
2379  // elements attached to them. Store that info. in the corresponding
2380  // data structures and be ready to send that info. to a root
2381  // processor. The root processor will be in charge of computing the
2382  // boundary coordinates for each segment of the boundary.
2383 
2384  // For each segment store the following information
2385  // --------------------------------------------------------------------
2386  // Stores the "rank" of the processor to the left of each segment,
2387  // zero if there is no processor to the left which states that the
2388  // segment is the first one on the boundary
2389  Vector<unsigned> left_processor_plus_one(nsegments);
2390 
2391  // Stores the "rank" of the processor to the right of each segment,
2392  // zero if there is no processor to the right which states that the
2393  // segment is the last one on the boundary
2394  Vector<unsigned> right_processor_plus_one(nsegments);
2395 
2396  // The id. of the halo element to the left of the segment, note that
2397  // this info. is not necessary if there is no processor to the left
2398  // of the segment
2399  Vector<unsigned> left_halo_element(nsegments);
2400 
2401  // The id. of the halo element to the right of the segment, note that
2402  // this info. is not necessary if there is no processor to the right
2403  // of the segment
2404  Vector<unsigned> right_halo_element(nsegments);
2405 
2406  // The id. of the haloed element to the left of the segment, note that
2407  // this info. is not necessary if there is no processor to the left
2408  // of the segment
2409  Vector<unsigned> left_haloed_element(nsegments);
2410 
2411  // The id. of the haloed element to the right of the segment, note
2412  // that this info. is not necessary if there is no processor to the
2413  // right of the segment
2414  Vector<unsigned> right_haloed_element(nsegments);
2415 
2416  // Go through all the segments and get the info.
2417  for (unsigned is = 0; is < nsegments; is++)
2418  {
2419  // Get access to the left most face element on the segment
2420  FiniteElement* left_face_ele_pt = segment_sorted_ele_pt[is].front();
2421 
2422  // Get the corresponding bulk element and check whether it is a halo
2423  // element or not
2424  FiniteElement* tmp_left_bulk_ele_pt =
2425  face_to_bulk_element_pt[left_face_ele_pt];
2426 
2427  // Check if the bulk element is halo
2428  if (tmp_left_bulk_ele_pt->is_halo())
2429  {
2430  // Then store the corresponding info.
2431  int left_proc = tmp_left_bulk_ele_pt->non_halo_proc_ID();
2432 #ifdef PARANOID
2433  if (left_proc < 0)
2434  {
2435  std::ostringstream error_message;
2436  error_message
2437  << "The current bulk element (left) is marked as halo but "
2438  << "the processor holding\nthe non-halo counterpart is "
2439  << "negative!\n";
2440  throw OomphLibError(error_message.str(),
2441  OOMPH_CURRENT_FUNCTION,
2442  OOMPH_EXCEPTION_LOCATION);
2443  }
2444 #endif
2445  // The processor "rank" to the left
2446  unsigned left_processor = static_cast<unsigned>(left_proc);
2447  left_processor_plus_one[is] = left_processor + 1;
2448 
2449  // Now get the id of the halo element to the left
2450  GeneralisedElement* left_element_pt = tmp_left_bulk_ele_pt;
2451 
2452  // Get the halo elements with left processor
2453  Vector<GeneralisedElement*> left_halo_element_pt =
2454  this->halo_element_pt(left_processor);
2455 
2456 #ifdef PARANOID
2457  // Flag to state that the halo element was found
2458  bool left_halo_element_found = false;
2459 #endif
2460 
2461  const unsigned n_halo_left = left_halo_element_pt.size();
2462  for (unsigned lh = 0; lh < n_halo_left; lh++)
2463  {
2464  if (left_element_pt == left_halo_element_pt[lh])
2465  {
2466  left_halo_element[is] = lh;
2467 #ifdef PARANOID
2468  left_halo_element_found = true;
2469 #endif
2470  break;
2471  }
2472  } // for (lh < n_halo_left)
2473 
2474 #ifdef PARANOID
2475  if (!left_halo_element_found)
2476  {
2477  std::ostringstream error_message;
2478  error_message
2479  << "The current bulk element (left) marked as halo was "
2480  << "not found in the vector of halo\nelements associated "
2481  << "with the (" << left_processor << ") processor.\n\n";
2482  throw OomphLibError(error_message.str(),
2483  OOMPH_CURRENT_FUNCTION,
2484  OOMPH_EXCEPTION_LOCATION);
2485  } // if (!left_halo_element_found)
2486 #endif
2487 
2488  // Get the left-most nonhalo element (use the backup list of
2489  // nonhalo elements)
2490  left_face_ele_pt = segment_sorted_nonhalo_ele_pt[is].front();
2491 
2492  // Get the corresponding bulk element
2493  tmp_left_bulk_ele_pt = face_to_bulk_element_pt[left_face_ele_pt];
2494 
2495 #ifdef PARANOID
2496  // This element should not be marked as halo
2497  if (tmp_left_bulk_ele_pt->is_halo())
2498  {
2499  std::ostringstream error_message;
2500  error_message
2501  << "The bulk element represetation of the left-most nonhalo face\n"
2502  << "element of the current segment (" << is
2503  << ") is marked as halo,\n"
2504  << "but the face element created from it is nonhalo\n";
2505  throw OomphLibError(error_message.str(),
2506  OOMPH_CURRENT_FUNCTION,
2507  OOMPH_EXCEPTION_LOCATION);
2508  } // if (tmp_left_bulk_ele_pt->is_halo())
2509 #endif
2510 
2511  // Cast from "FiniteElement*" to "GeneralisedElement*" to be able
2512  // to search in the haloed vector
2513  left_element_pt = tmp_left_bulk_ele_pt;
2514 
2515 #ifdef PARANOID
2516  // Flag to state that the haloed element was found
2517  bool left_haloed_element_found = false;
2518 #endif
2519 
2520  // Now get the id for the haloed element to the left, get the
2521  // haloed elements from the processor to the left
2522  Vector<GeneralisedElement*> left_haloed_element_pt =
2523  this->haloed_element_pt(left_processor);
2524 
2525  const unsigned nhaloed_left = left_haloed_element_pt.size();
2526  for (unsigned lhd = 0; lhd < nhaloed_left; lhd++)
2527  {
2528  if (left_element_pt == left_haloed_element_pt[lhd])
2529  {
2530  left_haloed_element[is] = lhd;
2531 #ifdef PARANOID
2532  left_haloed_element_found = true;
2533 #endif
2534  break;
2535  }
2536  } // for (lhd < nhaloed_left)
2537 
2538 #ifdef PARANOID
2539  if (!left_haloed_element_found)
2540  {
2541  std::ostringstream error_message;
2542  error_message
2543  << "The current bulk element (left) marked as haloed was "
2544  << "not found in the vector of haloed\nelements associated "
2545  << "with processor (" << left_processor << ").\n";
2546  throw OomphLibError(error_message.str(),
2547  OOMPH_CURRENT_FUNCTION,
2548  OOMPH_EXCEPTION_LOCATION);
2549  }
2550 #endif
2551  } // if (tmp_left_bulk_ele_pt->is_halo())
2552  else
2553  {
2554  // If not halo then state the info. to indicate that
2555  left_processor_plus_one[is] = 0;
2556  // Null this info.
2557  left_halo_element[is] = 0;
2558  // Null this info.
2559  left_haloed_element[is] = 0;
2560  }
2561 
2562  // Get access to the right most face element on the segment
2563  FiniteElement* right_face_ele_pt = segment_sorted_ele_pt[is].back();
2564 
2565  // Get the corresponding bulk element and check whether it is
2566  // a halo element or not
2567  FiniteElement* tmp_right_bulk_ele_pt =
2568  face_to_bulk_element_pt[right_face_ele_pt];
2569 
2570  // Check if the bulk element is halo
2571  if (tmp_right_bulk_ele_pt->is_halo())
2572  {
2573  // Then store the corresponding info.
2574  int right_proc = tmp_right_bulk_ele_pt->non_halo_proc_ID();
2575 #ifdef PARANOID
2576  if (right_proc < 0)
2577  {
2578  std::ostringstream error_message;
2579  error_message
2580  << "The current bulk element (right) is marked as halo but "
2581  << "the processor holding\nthe non-halo counterpart is "
2582  << "negative!\n";
2583  throw OomphLibError(error_message.str(),
2584  "TriangleMesh::compute_boundary_segments_"
2585  "connectivity_and_initial_zeta_values()",
2586  OOMPH_EXCEPTION_LOCATION);
2587  }
2588 #endif
2589  // The processor "rank" to the right
2590  unsigned right_processor = static_cast<unsigned>(right_proc);
2591  right_processor_plus_one[is] = right_processor + 1;
2592 
2593  // Now get the id of the halo element to the right
2594  GeneralisedElement* right_element_pt = tmp_right_bulk_ele_pt;
2595 
2596  // Get the halo elements with right processor
2597  Vector<GeneralisedElement*> right_halo_element_pt =
2598  this->halo_element_pt(right_processor);
2599 
2600 #ifdef PARANOID
2601  // Flag to state that the halo element was found
2602  bool right_halo_element_found = false;
2603 #endif
2604 
2605  const unsigned nhalo_right = right_halo_element_pt.size();
2606  for (unsigned rh = 0; rh < nhalo_right; rh++)
2607  {
2608  if (right_element_pt == right_halo_element_pt[rh])
2609  {
2610  right_halo_element[is] = rh;
2611 #ifdef PARANOID
2612  right_halo_element_found = true;
2613 #endif
2614  break;
2615  }
2616  } // for (rh < nhalo_right)
2617 #ifdef PARANOID
2618  if (!right_halo_element_found)
2619  {
2620  std::ostringstream error_message;
2621  error_message
2622  << "The current bulk element (right) marked as halo was not "
2623  << "found in the vector of halo\nelements associated with "
2624  << "the (" << right_processor << ") processor.\n\n";
2625  throw OomphLibError(error_message.str(),
2626  "TriangleMesh::compute_boundary_segments_"
2627  "connectivity_and_initial_zeta_values()",
2628  OOMPH_EXCEPTION_LOCATION);
2629  }
2630 #endif
2631 
2632  // Get the right-most nonhalo element (use the backup list of
2633  // nonhalo elements)
2634  right_face_ele_pt = segment_sorted_nonhalo_ele_pt[is].back();
2635 
2636  // Get the corresponding bulk element
2637  tmp_right_bulk_ele_pt = face_to_bulk_element_pt[right_face_ele_pt];
2638 #ifdef PARANOID
2639  // This element should not be marked as halo
2640  if (tmp_right_bulk_ele_pt->is_halo())
2641  {
2642  std::ostringstream error_message;
2643  error_message
2644  << "The bulk element represetation of the right-most nonhalo face\n"
2645  << "element of the current segment (" << is
2646  << ") is marked as halo,\n"
2647  << "but the face element created from it is nonhalo\n";
2648  throw OomphLibError(error_message.str(),
2649  "TriangleMesh::compute_boundary_segments_"
2650  "connectivity_and_initial_zeta_values()",
2651  OOMPH_EXCEPTION_LOCATION);
2652  } // if (tmp_right_bulk_ele_pt->is_halo())
2653 #endif
2654 
2655  // Cast from "FiniteElement*" to "GeneralisedElement*" to be able
2656  // to search in the haloed vector
2657  right_element_pt = tmp_right_bulk_ele_pt;
2658 
2659 #ifdef PARANOID
2660  // Flag to state that the haloed element was found
2661  bool right_haloed_element_found = false;
2662 #endif
2663 
2664  // Now get the id for the haloed element to the right
2665  Vector<GeneralisedElement*> right_haloed_element_pt =
2666  this->haloed_element_pt(right_processor);
2667 
2668  const unsigned nhaloed_right = right_haloed_element_pt.size();
2669  for (unsigned rhd = 0; rhd < nhaloed_right; rhd++)
2670  {
2671  if (right_element_pt == right_haloed_element_pt[rhd])
2672  {
2673  right_haloed_element[is] = rhd;
2674 #ifdef PARANOID
2675  right_haloed_element_found = true;
2676 #endif
2677  break;
2678  }
2679  } // for (rhd < nhaloed_right)
2680 
2681 #ifdef PARANOID
2682  if (!right_haloed_element_found)
2683  {
2684  std::ostringstream error_message;
2685  error_message
2686  << "The current bulk element (right) marked as haloed was not "
2687  << "found in the vector of haloed\nelements associated with "
2688  << "the (" << right_processor << ") processor.\n\n";
2689  throw OomphLibError(error_message.str(),
2690  "TriangleMesh::compute_boundary_segments_"
2691  "connectivity_and_initial_zeta_values()",
2692  OOMPH_EXCEPTION_LOCATION);
2693  }
2694 #endif
2695 
2696  } // if (tmp_right_bulk_ele_pt->is_halo())
2697  else
2698  {
2699  // If not halo then state the info. to indicate that
2700  right_processor_plus_one[is] = 0;
2701  // Null this info.
2702  right_halo_element[is] = 0;
2703  // Null this info.
2704  right_haloed_element[is] = 0;
2705  }
2706 
2707  } // for (is < nsegments). Used to get the halo info. of the
2708  // segments
2709 
2710  // Now we have all the info. to be sent to the root processor and
2711  // compute the correct (global) boundary coordinates for the current
2712  // boundary
2713 
2714  // The root processor will be in charge of performing the computing
2715  // of the coordinate values along the boundary, all the other
2716  // processors only send their info. and wait for receiving the new
2717  // starting values for each of its segments
2718 
2719  // Choose the root processor
2720  const unsigned root_processor = 0;
2721  // ------------------------------------------------------------------
2722  // Starts the MPI stage
2723 
2724  // The root processor receives the number of segments of each
2725  // processor associated to the current boundary
2726  Vector<unsigned> root_nsegments_per_processor(nproc);
2727  unsigned nsegments_mpi = nsegments;
2728  MPI_Gather(&nsegments_mpi,
2729  1,
2730  MPI_UNSIGNED,
2731  &root_nsegments_per_processor[0],
2732  1,
2733  MPI_UNSIGNED,
2734  root_processor,
2735  comm_pt->mpi_comm());
2736 
2737  // Package the info. and prepare it to be sent
2738  // For the packaged info. we send 7 data per each segment, the indexes
2739  // are as follow; 0 left proc, 1 right proc, 2 left halo, 3 right
2740  // halo, 4 left haloed, 5 right haloed and 6 for nvertices per
2741  // segment
2742  // The size of the package (unsigned)
2743  const unsigned spu = 7;
2744  Vector<unsigned> flat_packed_unsigned_send_data(nsegments * spu);
2745  for (unsigned is = 0; is < nsegments; is++)
2746  {
2747  flat_packed_unsigned_send_data[(spu * is) + 0] =
2748  left_processor_plus_one[is];
2749  flat_packed_unsigned_send_data[(spu * is) + 1] =
2750  right_processor_plus_one[is];
2751  flat_packed_unsigned_send_data[(spu * is) + 2] = left_halo_element[is];
2752  flat_packed_unsigned_send_data[(spu * is) + 3] = right_halo_element[is];
2753  flat_packed_unsigned_send_data[(spu * is) + 4] = left_haloed_element[is];
2754  flat_packed_unsigned_send_data[(spu * is) + 5] = right_haloed_element[is];
2755  flat_packed_unsigned_send_data[(spu * is) + 6] =
2756  nvertices_per_segment[is];
2757  }
2758 
2759  // How many data will this processor send
2760  const unsigned nudata_to_send = flat_packed_unsigned_send_data.size();
2761 
2762  // How many data does the root processor will receive from each
2763  // processor
2764  Vector<int> root_nudata_to_receive(nproc, 0);
2765  // Total number of data to receive from all processors
2766  unsigned root_nutotal_data_receive = 0;
2767  for (unsigned ip = 0; ip < nproc; ip++)
2768  {
2769  // Compute the number of data the root processor will receive from
2770  // each processor
2771  root_nudata_to_receive[ip] = root_nsegments_per_processor[ip] * spu;
2772  // Add on the total number of data to receive
2773  root_nutotal_data_receive += root_nudata_to_receive[ip];
2774  }
2775 
2776  // Stores and compute the offsets (in root) for the data received
2777  // from each processor
2778  Vector<int> root_uoffsets_receive(nproc, 0);
2779  root_uoffsets_receive[0] = 0;
2780  for (unsigned ip = 1; ip < nproc; ip++)
2781  {
2782  // Compute the offset to store the values from each processor
2783  root_uoffsets_receive[ip] =
2784  root_uoffsets_receive[ip - 1] + root_nudata_to_receive[ip - 1];
2785  }
2786 
2787  // Create at least one entry so we don't get a seg fault below
2788  if (flat_packed_unsigned_send_data.size() == 0)
2789  {
2790  flat_packed_unsigned_send_data.resize(1);
2791  }
2792 
2793  // Vector where to receive the info.
2794  Vector<unsigned> flat_packed_unsigned_receive_data(
2795  root_nutotal_data_receive);
2796  if (my_rank != root_processor)
2797  {
2798  // Create at least one entry so we don't get a seg fault below
2799  if (flat_packed_unsigned_receive_data.size() == 0)
2800  {
2801  flat_packed_unsigned_receive_data.resize(1);
2802  }
2803  } // if (my_rank!=root_processor)
2804 
2805  MPI_Gatherv(&flat_packed_unsigned_send_data[0], // Flat package to
2806  // send info. from
2807  // each processor
2808  nudata_to_send, // Total number of data to send from
2809  // each processor
2810  MPI_UNSIGNED,
2811  &flat_packed_unsigned_receive_data[0], // Container
2812  // where to
2813  // receive the
2814  // info. from all
2815  // the processors
2816  &root_nudata_to_receive[0], // Number of data to receive
2817  // from each processor
2818  &root_uoffsets_receive[0], // The offset to store the
2819  // info. from each processor
2820  MPI_UNSIGNED,
2821  root_processor, // The processor that receives all the
2822  // info.
2823  comm_pt->mpi_comm());
2824 
2825  // Clear the flat package to send
2826  flat_packed_unsigned_send_data.clear();
2827  flat_packed_unsigned_send_data.resize(0);
2828 
2829  // Package the info. and prepare it to be sent
2830  // For the packaged info. we send 1 data per each segment which is
2831  // at the moment the arclength of each segment
2832  // The size of the package
2833  const unsigned spd = 1;
2834  Vector<double> flat_packed_double_send_data(nsegments * spd);
2835  for (unsigned is = 0; is < nsegments; is++)
2836  {
2837  flat_packed_double_send_data[(spd * is) + 0] = segment_arclength[is];
2838  }
2839 
2840  // How many data will this processor send
2841  const unsigned nddata_to_send = flat_packed_double_send_data.size();
2842  // How many data does the root processor will receive from each
2843  // processor
2844  Vector<int> root_nddata_to_receive(nproc, 0);
2845  // Total number of data to receive from all processors
2846  unsigned root_ndtotal_data_receive = 0;
2847  for (unsigned ip = 0; ip < nproc; ip++)
2848  {
2849  root_nddata_to_receive[ip] = root_nsegments_per_processor[ip] * spd;
2850  root_ndtotal_data_receive += root_nddata_to_receive[ip];
2851  }
2852 
2853  // Stores and compute the offsets for the data received from each
2854  // processor
2855  Vector<int> root_doffsets_receive(nproc, 0);
2856  root_doffsets_receive[0] = 0;
2857  for (unsigned ip = 1; ip < nproc; ip++)
2858  {
2859  // Compute the offset to store the values from each processor
2860  root_doffsets_receive[ip] =
2861  root_doffsets_receive[ip - 1] + root_nddata_to_receive[ip - 1];
2862  }
2863 
2864  // Create at least one entry so we don't get a seg fault below
2865  if (flat_packed_double_send_data.size() == 0)
2866  {
2867  flat_packed_double_send_data.resize(1);
2868  }
2869 
2870  // Vector where to receive the info.
2871  Vector<double> flat_packed_double_receive_data(root_ndtotal_data_receive);
2872  if (my_rank != root_processor)
2873  {
2874  // Create at least one entry so we don't get a seg fault below
2875  if (flat_packed_double_receive_data.size() == 0)
2876  {
2877  flat_packed_double_receive_data.resize(1);
2878  }
2879  }
2880 
2881  MPI_Gatherv(&flat_packed_double_send_data[0], // Flat package to
2882  // send info. from
2883  // each processor
2884  nddata_to_send, // Total number of data to send from
2885  // each processor
2886  MPI_DOUBLE,
2887  &flat_packed_double_receive_data[0], // Container where
2888  // to receive the
2889  // info. from all
2890  // the processors
2891  &root_nddata_to_receive[0], // Number of data to receive
2892  // from each processor
2893  &root_doffsets_receive[0], // The offset to store the
2894  // info. from each processor
2895  MPI_DOUBLE,
2896  root_processor, // The processor that receives all the
2897  // info.
2898  comm_pt->mpi_comm());
2899 
2900  // Clear the flat package to send
2901  flat_packed_double_send_data.clear();
2902  flat_packed_double_send_data.resize(0);
2903 
2904  // The next three containers are only used by the root processor at
2905  // the end of its computations but it is necessary that all the
2906  // processors know them when calling back the info.
2907 
2908  // Container that state the initial arclength for each segments
2909  // of each processor
2910  Vector<Vector<double>> root_initial_segment_arclength(nproc);
2911 
2912  // Container that state the number of vertices before each segment
2913  // in a given processor
2914  Vector<Vector<unsigned>> root_nvertices_before_segment(nproc);
2915 
2916  // The root processor needs to tell the other processor if it was
2917  // necessary to reverse a segment. Each processor should therefore
2918  // invert the face elements that compose every segment that was
2919  // inverted by the root processor
2920  Vector<Vector<unsigned>> root_segment_inverted(nproc);
2921 
2922  // Used to store the accumulated arclength, used at the end of
2923  // communications to store the total arclength
2924  double root_accumulated_arclength = 0.0;
2925 
2926  // Store the accumulated number of vertices, it means the total number
2927  // of vertices before each segment (counter)
2928  unsigned root_accumulated_vertices_before_segment = 0;
2929 
2930  // The root processor is in charge of performing the connections
2931  // of the segments that define the complete boundary
2932  if (my_rank == root_processor)
2933  {
2934  // From the flat packaged received data re-create the data
2935  // structures storing the info. regarding the connectivity of the
2936  // segments, the number of vertices per segment and the local
2937  // arclength of each segment
2938 
2939  // Stores the "rank" of the processor to the left of each segment,
2940  // zero if there is no processor to the left which states that the
2941  // segment is the first one on the boundary
2942  Vector<Vector<unsigned>> root_left_processor_plus_one(nproc);
2943 
2944  // Stores the "rank" of the processor to the right of each segment,
2945  // zero if there is no processor to the right which states that the
2946  // segment is the last one on the boundary
2947  Vector<Vector<unsigned>> root_right_processor_plus_one(nproc);
2948 
2949  // The id. of the halo element to the left of the segment, note that
2950  // this info. is not necessary if there is no processor to the left
2951  // of the segment or if the processor has no info about the boundary
2952  Vector<Vector<unsigned>> root_left_halo_element(nproc);
2953 
2954  // The id. of the halo element to the right of the segment, note
2955  // that this info. is not necessary if there is no processor to
2956  // the right of the segment or if the processor has no info about
2957  // the boundary
2958  Vector<Vector<unsigned>> root_right_halo_element(nproc);
2959 
2960  // The id. of the haloed element to the left of the segment, note
2961  // that this info. is not necessary if there is no processor to
2962  // the left of the segment or if the processor has no info about
2963  // the boundary
2964  Vector<Vector<unsigned>> root_left_haloed_element(nproc);
2965 
2966  // The id. of the haloed element to the right of the segment, note
2967  // that this info. is not necessary if there is no processor to the
2968  // right of the segment or if the processor has no info about the
2969  // boundary
2970  Vector<Vector<unsigned>> root_right_haloed_element(nproc);
2971 
2972  // The number of vertices per segment in each processor
2973  Vector<Vector<unsigned>> root_nvertices_per_segment(nproc);
2974 
2975  // The arclength of each of the segments in the processors
2976  Vector<Vector<double>> root_segment_arclength(nproc);
2977 
2978  unsigned ucounter = 0;
2979  unsigned dcounter = 0;
2980  for (unsigned ip = 0; ip < nproc; ip++)
2981  {
2982  // Get the number of segments in the current processor
2983  const unsigned nsegs_iproc = root_nsegments_per_processor[ip];
2984 
2985  root_left_processor_plus_one[ip].resize(nsegs_iproc);
2986  root_right_processor_plus_one[ip].resize(nsegs_iproc);
2987  root_left_halo_element[ip].resize(nsegs_iproc);
2988  root_right_halo_element[ip].resize(nsegs_iproc);
2989  root_left_haloed_element[ip].resize(nsegs_iproc);
2990  root_right_haloed_element[ip].resize(nsegs_iproc);
2991 
2992  // Additional info.
2993  root_nvertices_per_segment[ip].resize(nsegs_iproc);
2994  root_segment_arclength[ip].resize(nsegs_iproc);
2995  root_segment_inverted[ip].resize(nsegs_iproc);
2996 
2997  // Extract the info. from the BIG package received from all
2998  // processors
2999  for (unsigned is = 0; is < nsegs_iproc; is++)
3000  {
3001  // ------ The flat unsigned package ------
3002  root_left_processor_plus_one[ip][is] =
3003  flat_packed_unsigned_receive_data[ucounter++];
3004  root_right_processor_plus_one[ip][is] =
3005  flat_packed_unsigned_receive_data[ucounter++];
3006  root_left_halo_element[ip][is] =
3007  flat_packed_unsigned_receive_data[ucounter++];
3008  root_right_halo_element[ip][is] =
3009  flat_packed_unsigned_receive_data[ucounter++];
3010  root_left_haloed_element[ip][is] =
3011  flat_packed_unsigned_receive_data[ucounter++];
3012  root_right_haloed_element[ip][is] =
3013  flat_packed_unsigned_receive_data[ucounter++];
3014  root_nvertices_per_segment[ip][is] =
3015  flat_packed_unsigned_receive_data[ucounter++];
3016 
3017  // ------ The flat double package ------
3018  root_segment_arclength[ip][is] =
3019  flat_packed_double_receive_data[dcounter++];
3020  } // for (is < nsegs_iproc)
3021  } // for (ip < nproc)
3022 
3023  // Now the root processor has all the info. to find out the
3024  // CONNECTIVITY of the segments in each processor
3025 
3026  // Container that stores the info. related with the connectivity
3027  // of the segments of each processor
3028  Vector<Vector<int>> left_connected_segment_plus_one(nproc);
3029  Vector<Vector<int>> right_connected_segment_plus_one(nproc);
3030  for (unsigned ip = 0; ip < nproc; ip++)
3031  {
3032  const unsigned nsegs_iproc = root_nsegments_per_processor[ip];
3033  left_connected_segment_plus_one[ip].resize(nsegs_iproc, -1);
3034  right_connected_segment_plus_one[ip].resize(nsegs_iproc, -1);
3035  } // for (ip < nprocs)
3036 
3037  // In charge of storing the connectivity of the segments, the pair
3038  // indicates the processor and the segment number
3039  std::list<std::pair<unsigned, unsigned>> proc_seg_connectivity;
3040  proc_seg_connectivity.clear();
3041 
3042  // Done segments on processor
3043  std::map<std::pair<unsigned, unsigned>, bool> done_segment;
3044 
3045  // Take the first segment of the first processor with segments and
3046  // add it to the list of segments
3047  unsigned left_proc = 0;
3048  unsigned right_proc = 0;
3049  unsigned left_seg = 0;
3050  unsigned right_seg = 0;
3051  for (unsigned ip = 0; ip < nproc; ip++)
3052  {
3053  const unsigned nsegs_iproc = root_nsegments_per_processor[ip];
3054  if (nsegs_iproc > 0)
3055  {
3056  right_proc = left_proc = ip;
3057  right_seg = left_seg = 0;
3058  break; // Break because it is the first processor with at
3059  // least one segment
3060  }
3061  } // for (ip < nproc)
3062 
3063  // ... and add it to the list of segments
3064  std::pair<unsigned, unsigned> add_segment =
3065  std::make_pair(left_proc, left_seg);
3066  done_segment[add_segment] = true;
3067  proc_seg_connectivity.push_back(add_segment);
3068 
3069  // Flags to indicate when a segment was added to the left or right
3070  // of the current list of segments
3071  bool added_segment_to_the_left = false;
3072  bool added_segment_to_the_right = false;
3073 
3074  do // while(added_segment_to_the_left || added_segment_to_the_right)
3075  {
3076  // Read the left-most processor and segment in the list
3077  std::pair<unsigned, unsigned> left_pair = proc_seg_connectivity.front();
3078  left_proc = left_pair.first;
3079  left_seg = left_pair.second;
3080 
3081  // Get the processor number to the left of the left-most
3082  // segment in the list
3083  const unsigned new_left_proc =
3084  root_left_processor_plus_one[left_proc][left_seg];
3085 
3086  if (new_left_proc != 0)
3087  {
3088  // Initialise flag
3089  added_segment_to_the_left = false;
3090  // Get the left halo element id
3091  const unsigned left_halo_id =
3092  root_left_halo_element[left_proc][left_seg];
3093 
3094  // Get the left haloed element id
3095  const unsigned left_haloed_id =
3096  root_left_haloed_element[left_proc][left_seg];
3097 
3098  // Go through the segments on the new left processor and look
3099  // for the corresponding left_halo_id in the haloed_ids
3100  const unsigned nsegs_new_left_proc =
3101  root_nsegments_per_processor[new_left_proc - 1];
3102 
3103  for (unsigned ils = 0; ils < nsegs_new_left_proc; ils++)
3104  {
3105  std::pair<unsigned, unsigned> candidate_seg =
3106  std::make_pair(new_left_proc - 1, ils);
3107 
3108  // Check that the segment has not been already added
3109  if (!done_segment[candidate_seg])
3110  {
3111  // Only consider the segments on new left processor which
3112  // right processor is the current one (left_proc)
3113  const unsigned right_proc_of_new_left_proc =
3114  root_right_processor_plus_one[new_left_proc - 1][ils];
3115  // Also get the left_proc_of_new_left_proc (in case that it
3116  // be necessary to invert the segment)
3117  const unsigned left_proc_of_new_left_proc =
3118  root_left_processor_plus_one[new_left_proc - 1][ils];
3119  // Check the not inverted case (to the left and not
3120  // inverted)
3121  if (right_proc_of_new_left_proc != 0 &&
3122  right_proc_of_new_left_proc - 1 == left_proc)
3123  {
3124  // Get the haloed/haloed element id of the current segment
3125  // in the new left processor and compare it to the
3126  // halo/haloed element id of the left_processor
3127  const unsigned right_halo_id =
3128  root_right_halo_element[new_left_proc - 1][ils];
3129  const unsigned right_haloed_id =
3130  root_right_haloed_element[new_left_proc - 1][ils];
3131  if (left_halo_id == right_haloed_id &&
3132  left_haloed_id == right_halo_id)
3133  {
3134  // We have a match of the segments (store the segment
3135  // number plus one on the processor to the left)
3136  left_connected_segment_plus_one[left_proc][left_seg] =
3137  ils + 1;
3138  // Add the pair to the connectivity list
3139  proc_seg_connectivity.push_front(candidate_seg);
3140  added_segment_to_the_left = true;
3141  break;
3142  }
3143  } // if (right_proc_of_new_left_proc-1 == left_proc)
3144 
3145  // Check the inverted case (to the left and inverted)
3146  if (left_proc_of_new_left_proc != 0 &&
3147  left_proc_of_new_left_proc - 1 == left_proc)
3148  {
3149  // Get the haloed element id of the current segment
3150  // (inverted version) in the new left processor and
3151  // compare it to the halo element id of the left_processor
3152  const unsigned inv_left_halo_id =
3153  root_left_halo_element[new_left_proc - 1][ils];
3154  const unsigned inv_left_haloed_id =
3155  root_left_haloed_element[new_left_proc - 1][ils];
3156  if (left_halo_id == inv_left_haloed_id &&
3157  left_haloed_id == inv_left_halo_id)
3158  {
3159  // We have a match of the segments (store the segment
3160  // number plus one on the processor to the left)
3161  left_connected_segment_plus_one[left_proc][left_seg] =
3162  ils + 1;
3163  // Add the pair to the connectivity list
3164  proc_seg_connectivity.push_front(candidate_seg);
3165 
3166  // In addition to the connectivity we need to invert the
3167  // segment (the information)
3168  const unsigned tmp_proc =
3169  root_left_processor_plus_one[new_left_proc - 1][ils];
3170  const unsigned tmp_halo =
3171  root_left_halo_element[new_left_proc - 1][ils];
3172  const unsigned tmp_haloed =
3173  root_left_haloed_element[new_left_proc - 1][ils];
3174 
3175  root_left_processor_plus_one[new_left_proc - 1][ils] =
3176  root_right_processor_plus_one[new_left_proc - 1][ils];
3177  root_left_halo_element[new_left_proc - 1][ils] =
3178  root_right_halo_element[new_left_proc - 1][ils];
3179  root_left_haloed_element[new_left_proc - 1][ils] =
3180  root_right_haloed_element[new_left_proc - 1][ils];
3181 
3182  root_right_processor_plus_one[new_left_proc - 1][ils] =
3183  tmp_proc;
3184  root_right_halo_element[new_left_proc - 1][ils] = tmp_halo;
3185  root_right_haloed_element[new_left_proc - 1][ils] =
3186  tmp_haloed;
3187 
3188  // ... and mark the segment as inverted in the root
3189  // processor to inform back to the owner processor
3190  root_segment_inverted[new_left_proc - 1][ils] = 1;
3191 
3192  added_segment_to_the_left = true;
3193  break;
3194  }
3195  } // if (left_proc_of_new_left_proc-1 == left_proc)
3196  } // if (!done_segment[candidate_segment])
3197  } // for (ils < nsegs_new_left_proc)
3198 
3199 #ifdef PARANOID
3200  if (!added_segment_to_the_left)
3201  {
3202  std::ostringstream error_message;
3203  error_message
3204  << "The corresponding processor and segment to the left of "
3205  << "the current left\nmost segment was not found\n";
3206  throw OomphLibError(error_message.str(),
3207  "TriangleMesh::compute_boundary_segments_"
3208  "connectivity_and_initial_zeta_values()",
3209  OOMPH_EXCEPTION_LOCATION);
3210  }
3211 #endif
3212  } // if (new_left_proc != 0)
3213  else
3214  {
3215  // No more segments to the left
3216  added_segment_to_the_left = false;
3217  }
3218 
3219  // Read the info. of the right processor and the right segment
3220  std::pair<unsigned, unsigned> right_pair = proc_seg_connectivity.back();
3221  right_proc = right_pair.first;
3222  right_seg = right_pair.second;
3223 
3224  // Get the processor number to the right of the right-most
3225  // segment in the list
3226  const unsigned new_right_proc =
3227  root_right_processor_plus_one[right_proc][right_seg];
3228 
3229  if (new_right_proc != 0)
3230  {
3231  // Initialise flag
3232  added_segment_to_the_right = false;
3233  // Get the right halo element id
3234  const unsigned right_halo_id =
3235  root_right_halo_element[right_proc][right_seg];
3236 
3237  // Get the right halo element id
3238  const unsigned right_haloed_id =
3239  root_right_haloed_element[right_proc][right_seg];
3240 
3241  // Go through the segments on the new right processor and look
3242  // for the corresponding right_halo_id in the haloed_ids
3243  const unsigned nsegs_new_right_proc =
3244  root_nsegments_per_processor[new_right_proc - 1];
3245 
3246  for (unsigned irs = 0; irs < nsegs_new_right_proc; irs++)
3247  {
3248  std::pair<unsigned, unsigned> candidate_seg =
3249  std::make_pair(new_right_proc - 1, irs);
3250 
3251  // Check that the segment has not been already added
3252  if (!done_segment[candidate_seg])
3253  {
3254  // Only consider the segments on new right processor which
3255  // left processor is the current one (right_proc)
3256  const unsigned left_proc_of_new_right_proc =
3257  root_left_processor_plus_one[new_right_proc - 1][irs];
3258  // Also get the right_proc_of_new_right_proc (in case
3259  // that it be necessary to invert the segment)
3260  const unsigned right_proc_of_new_right_proc =
3261  root_right_processor_plus_one[new_right_proc - 1][irs];
3262  // Check the not inverted case (to the right and not
3263  // inverted)
3264  if (left_proc_of_new_right_proc != 0 &&
3265  left_proc_of_new_right_proc - 1 == right_proc)
3266  {
3267  // Get the haloed element id of the current segment in the
3268  // new right processor and compare it to the halo element
3269  // id of the right_processor
3270  const unsigned left_halo_id =
3271  root_left_halo_element[new_right_proc - 1][irs];
3272  const unsigned left_haloed_id =
3273  root_left_haloed_element[new_right_proc - 1][irs];
3274 
3275  if (right_halo_id == left_haloed_id &&
3276  right_haloed_id == left_halo_id)
3277  {
3278  // We have a match of the segments (store the segment
3279  // number plus one on the processor to the right)
3280  right_connected_segment_plus_one[right_proc][right_seg] =
3281  irs + 1;
3282  // Add the connectivity information to the list
3283  proc_seg_connectivity.push_back(candidate_seg);
3284  added_segment_to_the_right = true;
3285  break;
3286  }
3287  } // if (left_proc_of_new_right_proc-1 == right_proc)
3288 
3289  // Check the inverted case (to the right and inverted)
3290  if (right_proc_of_new_right_proc != 0 &&
3291  right_proc_of_new_right_proc - 1 == right_proc)
3292  {
3293  // Get the haloed element id of the current segment
3294  // (inverted version) in the new right processor and
3295  // compare it to the halo element id of the
3296  // right_processor
3297  const unsigned inv_right_halo_id =
3298  root_right_halo_element[new_right_proc - 1][irs];
3299  const unsigned inv_right_haloed_id =
3300  root_right_haloed_element[new_right_proc - 1][irs];
3301  if (right_halo_id == inv_right_haloed_id &&
3302  right_haloed_id == inv_right_halo_id)
3303  {
3304  // We have a match of the segments (store the segment
3305  // number plus one on the processor to the right)
3306  right_connected_segment_plus_one[right_proc][right_seg] =
3307  irs + 1;
3308  // Add the connectivity information to the list
3309  proc_seg_connectivity.push_back(candidate_seg);
3310  // In addition to the connectivity we need to invert the
3311  // segment
3312  const unsigned tmp_proc =
3313  root_left_processor_plus_one[new_right_proc - 1][irs];
3314  const unsigned tmp_halo =
3315  root_left_halo_element[new_right_proc - 1][irs];
3316  const unsigned tmp_haloed =
3317  root_left_haloed_element[new_right_proc - 1][irs];
3318 
3319  root_left_processor_plus_one[new_right_proc - 1][irs] =
3320  root_right_processor_plus_one[new_right_proc - 1][irs];
3321  root_left_halo_element[new_right_proc - 1][irs] =
3322  root_right_halo_element[new_right_proc - 1][irs];
3323  root_left_haloed_element[new_right_proc - 1][irs] =
3324  root_right_haloed_element[new_right_proc - 1][irs];
3325 
3326  root_right_processor_plus_one[new_right_proc - 1][irs] =
3327  tmp_proc;
3328  root_right_halo_element[new_right_proc - 1][irs] = tmp_halo;
3329  root_right_haloed_element[new_right_proc - 1][irs] =
3330  tmp_haloed;
3331 
3332  // ... and mark the segment as inverted in the root
3333  // processor to inform back to the owner processor
3334  root_segment_inverted[new_right_proc - 1][irs] = 1;
3335 
3336  added_segment_to_the_right = true;
3337  break;
3338  }
3339  } // if (right_proc_of_new_right_proc-1 == right_proc)
3340  } // if (!done_segment[candidate_segment])
3341  } // for (irs < nsegs_new_left_proc)
3342 
3343 #ifdef PARANOID
3344  if (!added_segment_to_the_right)
3345  {
3346  std::ostringstream error_message;
3347  error_message
3348  << "The corresponding processor and segment to the right of "
3349  << "the current right\nmost segment was not found\n";
3350  throw OomphLibError(error_message.str(),
3351  "TriangleMesh::compute_boundary_segments_"
3352  "connectivity_and_initial_zeta_values()",
3353  OOMPH_EXCEPTION_LOCATION);
3354  }
3355 #endif
3356  } // if (new_right_proc != 0)
3357  else
3358  {
3359  // No more segments to the left
3360  added_segment_to_the_right = false;
3361  }
3362 
3363  } while (added_segment_to_the_left || added_segment_to_the_right);
3364 
3365  // Once we have connected the segments then we can compute the
3366  // initial and final zeta values based on the arclength of each
3367  // individual segment
3368 
3369  // Get the total number of segments, which MUST be the same as the
3370  // total number of segments in all processors
3371  const unsigned ntotal_segments = proc_seg_connectivity.size();
3372 #ifdef PARANOID
3373  unsigned tmp_total_segments = 0;
3374  for (unsigned ip = 0; ip < nproc; ip++)
3375  {
3376  tmp_total_segments += root_nsegments_per_processor[ip];
3377  }
3378 
3379  // Check that the total number of segments in all processors is
3380  // the same as the number of segments that form the boundary
3381  if (ntotal_segments != tmp_total_segments)
3382  {
3383  std::ostringstream error_message;
3384  error_message << "The number of sorted segments (" << ntotal_segments
3385  << ") on "
3386  << "boundary (" << b
3387  << ")\nis different from the total number of "
3388  << "segments (" << tmp_total_segments
3389  << ") in all\nprocessors.\n\n";
3390  throw OomphLibError(error_message.str(),
3391  OOMPH_CURRENT_FUNCTION,
3392  OOMPH_EXCEPTION_LOCATION);
3393  } // if (ntotal_segments!=tmp_total_segments)
3394 #endif
3395 
3396  // Now that we know the connectivity of the segments we can
3397  // compute the initial arclength of each segment in the
3398  // processors. Additionally we also get the number of vertices
3399  // before each of the segments. Resize the containers considering
3400  // the number of segments in each processor
3401  for (unsigned ip = 0; ip < nproc; ip++)
3402  {
3403  const unsigned nsegs_iproc = root_nsegments_per_processor[ip];
3404  root_initial_segment_arclength[ip].resize(nsegs_iproc);
3405  root_nvertices_before_segment[ip].resize(nsegs_iproc);
3406  }
3407 
3408  Vector<double> aux_initial_segment_arclength(ntotal_segments);
3409  Vector<unsigned> aux_nvertices_before_segment(ntotal_segments);
3410 
3411  ucounter = 0;
3412  for (std::list<std::pair<unsigned, unsigned>>::iterator it_list =
3413  proc_seg_connectivity.begin();
3414  it_list != proc_seg_connectivity.end();
3415  it_list++)
3416  {
3417  const unsigned iproc = static_cast<unsigned>((*it_list).first);
3418  const unsigned iseg = static_cast<unsigned>((*it_list).second);
3419  const double iseg_arclength = root_segment_arclength[iproc][iseg];
3420  const unsigned iseg_nvertices = root_nvertices_per_segment[iproc][iseg];
3421 
3422  aux_initial_segment_arclength[ucounter] = root_accumulated_arclength;
3423  aux_nvertices_before_segment[ucounter] =
3424  root_accumulated_vertices_before_segment;
3425 
3426  // Set the initial zeta value for the segment
3427  root_initial_segment_arclength[iproc][iseg] =
3428  root_accumulated_arclength;
3429  // Set the number of vertices before the current segment
3430  root_nvertices_before_segment[iproc][iseg] =
3431  root_accumulated_vertices_before_segment;
3432 
3433  // Add the arclength of the segment to the global arclength
3434  root_accumulated_arclength += iseg_arclength;
3435  // Add the number of vertices to the global number of vertices
3436  root_accumulated_vertices_before_segment += iseg_nvertices - 1;
3437 
3438  // Increase the counter
3439  ucounter++;
3440  } // for (loop over the sorted segments to assigne initial
3441  // arlength and initial number of vertices)
3442 
3443  // Increase by one to get the total number of vertices on the
3444  // boundary
3445  root_accumulated_vertices_before_segment++;
3446 
3447  // Get the processors with the initial and final segment.
3448  proc_with_initial_seg = proc_seg_connectivity.front().first;
3449  proc_with_final_seg = proc_seg_connectivity.back().first;
3450  // Also get the corresponding initial and final segment indexes
3451  // (on the initial and final processors)
3452  initial_segment = proc_seg_connectivity.front().second;
3453  final_segment = proc_seg_connectivity.back().second;
3454 
3455  } // if (my_rank == root_processor)
3456 
3457  // Get the total number of segments
3458  unsigned root_ntotal_segments = 0;
3459  for (unsigned ip = 0; ip < nproc; ip++)
3460  {
3461  root_ntotal_segments += root_nsegments_per_processor[ip];
3462  }
3463 
3464  // Package the info. that will be sent to each processor. For the
3465  // unsigned package we send the number of vertices before each
3466  // segment in each processor and whether it was inverted or not
3467  // Package size
3468  const unsigned rspu = 2;
3469  flat_packed_unsigned_send_data.clear();
3470  flat_packed_unsigned_send_data.resize(root_ntotal_segments * rspu);
3471  unsigned ucounter = 0;
3472  // Collect the info. from all the segments in the processors
3473  for (unsigned ip = 0; ip < nproc; ip++)
3474  {
3475  const unsigned nsegs_iproc = root_nsegments_per_processor[ip];
3476  for (unsigned is = 0; is < nsegs_iproc; is++)
3477  {
3478  flat_packed_unsigned_send_data[ucounter++] =
3479  root_nvertices_before_segment[ip][is];
3480  flat_packed_unsigned_send_data[ucounter++] =
3481  root_segment_inverted[ip][is];
3482  } // for (is < nsegs_iproc)
3483  } // for (ip < nproc)
3484 
3485  // How many data does the root processor will send to each processor
3486  Vector<int> root_nudata_to_send(nproc, 0);
3487  for (unsigned ip = 0; ip < nproc; ip++)
3488  {
3489  // Get the number of data to send to ip processor
3490  root_nudata_to_send[ip] = root_nsegments_per_processor[ip] * rspu;
3491  }
3492 
3493  // Store and compute the offsets for the data sent to each processor
3494  Vector<int> root_uoffsets_send(nproc, 0);
3495  root_uoffsets_send[0] = 0;
3496  for (unsigned ip = 1; ip < nproc; ip++)
3497  {
3498  // Compute the offset to send the values to each processor
3499  root_uoffsets_send[ip] =
3500  root_uoffsets_send[ip - 1] + root_nudata_to_send[ip - 1];
3501  }
3502 
3503  // Number of data to receive from root
3504  unsigned nutotal_data_receive = nsegments * rspu;
3505 
3506  if (my_rank != root_processor)
3507  {
3508  // Create at least one entry so we don't get a seg fault below
3509  if (flat_packed_unsigned_send_data.size() == 0)
3510  {
3511  flat_packed_unsigned_send_data.resize(1);
3512  }
3513  }
3514 
3515  // Clear and resize the vector where to receive the info.
3516  flat_packed_unsigned_receive_data.clear();
3517  flat_packed_unsigned_receive_data.resize(nutotal_data_receive);
3518  // Create at least one entry so we don't get a seg fault below
3519  if (flat_packed_unsigned_receive_data.size() == 0)
3520  {
3521  flat_packed_unsigned_receive_data.resize(1);
3522  }
3523 
3524  MPI_Scatterv(&flat_packed_unsigned_send_data[0],
3525  &root_nudata_to_send[0],
3526  &root_uoffsets_send[0],
3527  MPI_UNSIGNED,
3528  &flat_packed_unsigned_receive_data[0],
3529  nutotal_data_receive,
3530  MPI_UNSIGNED,
3531  root_processor,
3532  comm_pt->mpi_comm());
3533 
3534  // Package the info. that will be sent to each processor, for the
3535  // double package we send (one data per segment) the initial
3536  // arclength for each segment
3537  const unsigned rspd = 1;
3538  flat_packed_double_send_data.clear();
3539  flat_packed_double_send_data.resize(root_ntotal_segments * rspd);
3540  unsigned dcounter = 0;
3541  // Collect the info. from all the segments in the processors
3542  for (unsigned ip = 0; ip < nproc; ip++)
3543  {
3544  const unsigned nsegs_iproc = root_nsegments_per_processor[ip];
3545  for (unsigned is = 0; is < nsegs_iproc; is++)
3546  {
3547  flat_packed_double_send_data[dcounter++] =
3548  root_initial_segment_arclength[ip][is];
3549  }
3550  }
3551 
3552  // How many data does the root processor will send to each processor
3553  Vector<int> root_nddata_to_send(nproc, 0);
3554  for (unsigned ip = 0; ip < nproc; ip++)
3555  {
3556  // Number of data send to ip processor
3557  root_nddata_to_send[ip] = root_nsegments_per_processor[ip] * rspd;
3558  }
3559 
3560  // Store and compute the offsets for the data sent to each processor
3561  Vector<int> root_doffsets_send(nproc, 0);
3562  root_doffsets_send[0] = 0;
3563  for (unsigned ip = 1; ip < nproc; ip++)
3564  {
3565  // Compute the offset to send the values to each processor
3566  root_doffsets_send[ip] =
3567  root_doffsets_send[ip - 1] + root_nddata_to_send[ip - 1];
3568  }
3569 
3570  // Number of double data to receive from root
3571  unsigned ndtotal_data_receive = nsegments * rspd;
3572 
3573  if (my_rank != root_processor)
3574  {
3575  // Create at least one entry so we don't get a seg fault below
3576  if (flat_packed_double_send_data.size() == 0)
3577  {
3578  flat_packed_double_send_data.resize(1);
3579  }
3580  }
3581 
3582  // Clear and resize the vector where to receive the info.
3583  flat_packed_double_receive_data.clear();
3584  flat_packed_double_receive_data.resize(ndtotal_data_receive);
3585  // Create at least one entry so we don't get a seg fault below
3586  if (flat_packed_double_receive_data.size() == 0)
3587  {
3588  flat_packed_double_receive_data.resize(1);
3589  }
3590 
3591  MPI_Scatterv(&flat_packed_double_send_data[0],
3592  &root_nddata_to_send[0],
3593  &root_doffsets_send[0],
3594  MPI_DOUBLE,
3595  &flat_packed_double_receive_data[0],
3596  ndtotal_data_receive,
3597  MPI_DOUBLE,
3598  root_processor,
3599  comm_pt->mpi_comm());
3600 
3601  // Read if the segments need to be inverted and read the initial
3602  // arclengths
3603  ucounter = 0;
3604  dcounter = 0;
3605 
3606  // Read the info. from the flat package and store it in their
3607  // corresponding containers
3608  for (unsigned is = 0; is < nsegments; is++)
3609  {
3610  // The flat unsigned package
3611  nvertices_before_segment[is] =
3612  flat_packed_unsigned_receive_data[ucounter++];
3613  // The segment inverted flag
3614  segment_inverted[is] = flat_packed_unsigned_receive_data[ucounter++];
3615  // The flat double package
3616  initial_segment_arclength[is] =
3617  flat_packed_double_receive_data[dcounter++];
3618  } // for (is < nsegments)
3619 
3620  // Perform two additional communications to get the total number of
3621  // vertices, the processors with the initial and final segments, the
3622  // corresponding initial and final segments ...
3623  const unsigned numore_info = 5;
3624  Vector<unsigned> flat_package_unsigned_more_info(numore_info);
3625  // Prepare the info ...
3626  flat_package_unsigned_more_info[0] =
3627  root_accumulated_vertices_before_segment;
3628  flat_package_unsigned_more_info[1] = proc_with_initial_seg;
3629  flat_package_unsigned_more_info[2] = proc_with_final_seg;
3630  flat_package_unsigned_more_info[3] = initial_segment;
3631  flat_package_unsigned_more_info[4] = final_segment;
3632 
3633  // Send the info. to all processors
3634  MPI_Bcast(&flat_package_unsigned_more_info[0],
3635  numore_info,
3636  MPI_UNSIGNED,
3637  root_processor,
3638  comm_pt->mpi_comm());
3639 
3640  // ... and store the info. in the proper containers
3641  root_accumulated_vertices_before_segment =
3642  flat_package_unsigned_more_info[0];
3643  proc_with_initial_seg = flat_package_unsigned_more_info[1];
3644  proc_with_final_seg = flat_package_unsigned_more_info[2];
3645  initial_segment = flat_package_unsigned_more_info[3];
3646  final_segment = flat_package_unsigned_more_info[4];
3647 
3648  // Do the same for the maximum zeta value
3649  MPI_Bcast(&root_accumulated_arclength,
3650  1,
3651  MPI_DOUBLE,
3652  root_processor,
3653  comm_pt->mpi_comm());
3654 
3655  // -----------------------------------------------------------------
3656  // Clear the storage to store the data that will be used by the
3657  // setup boundary coordinates method, if we do not perform the
3658  // cleaning then previous data from previous iterations will remain
3659  // there
3660  // -----------------------------------------------------------------
3661  // The info. for the boundary
3662  Boundary_initial_coordinate[b].clear();
3663  Boundary_final_coordinate[b].clear();
3664 
3665  Boundary_initial_zeta_coordinate[b].clear();
3666  Boundary_final_zeta_coordinate[b].clear();
3667 
3668  // The info. for the segments
3669  Boundary_segment_inverted[b].clear();
3670  Boundary_segment_initial_coordinate[b].clear();
3671  Boundary_segment_final_coordinate[b].clear();
3672 
3673  Boundary_segment_initial_zeta[b].clear();
3674  Boundary_segment_final_zeta[b].clear();
3675 
3676  Boundary_segment_initial_arclength[b].clear();
3677  Boundary_segment_final_arclength[b].clear();
3678 
3679  // Now copy all the info. to the containers to be sent to any other
3680  // mesh (in the adaptation method)
3681  for (unsigned is = 0; is < nsegments; is++)
3682  {
3683  // At this point we can get the initial and final coordinates for
3684  // each segment
3685  Vector<double> first_seg_coord(2);
3686  Vector<double> last_seg_coord(2);
3687 
3688  // In order to get the first and last coordinates of each segment we
3689  // first need to identify the first and last nonhalo element of each
3690  // segment, and then get the first and last node of the segment
3691 
3692  // Get the first nonhalo face element on the segment
3693  FiniteElement* first_seg_ele_pt =
3694  segment_sorted_nonhalo_ele_pt[is].front();
3695 
3696 #ifdef PARANOID
3697  // Check if the face element is nonhalo, it shouldn't, but better
3698  // check
3699  if (first_seg_ele_pt->is_halo())
3700  {
3701  std::ostringstream error_message;
3702  error_message << "The first face element in the (" << is
3703  << ")-th segment is halo\n";
3704  throw OomphLibError(error_message.str(),
3705  "TriangleMesh::compute_boundary_segments_"
3706  "connectivity_and_initial_zeta_values()",
3707  OOMPH_EXCEPTION_LOCATION);
3708  } // if (tmp_first_bulk_ele_pt->is_halo())
3709 #endif
3710 
3711  // Number of nodes
3712  const unsigned nnod = first_seg_ele_pt->nnode();
3713 
3714  // Get the first node of the current segment
3715  Node* first_seg_node_pt = first_seg_ele_pt->node_pt(0);
3716  if (is_inverted[first_seg_ele_pt])
3717  {
3718  first_seg_node_pt = first_seg_ele_pt->node_pt(nnod - 1);
3719  }
3720 
3721  // Get the last nonhalo face element on the segment
3722  FiniteElement* last_seg_ele_pt = segment_sorted_nonhalo_ele_pt[is].back();
3723 
3724 #ifdef PARANOID
3725  // Check if the face element is nonhalo, it shouldn't, but better
3726  // check
3727  if (last_seg_ele_pt->is_halo())
3728  {
3729  std::ostringstream error_message;
3730  error_message << "The last face element in the (" << is
3731  << ")-th segment is halo\n";
3732  throw OomphLibError(error_message.str(),
3733  "TriangleMesh::compute_boundary_segments_"
3734  "connectivity_and_initial_zeta_values()",
3735  OOMPH_EXCEPTION_LOCATION);
3736  } // if (tmp_first_bulk_ele_pt->is_halo())
3737 #endif
3738 
3739  // Get the last node of the current segment
3740  Node* last_seg_node_pt = last_seg_ele_pt->node_pt(nnod - 1);
3741  if (is_inverted[last_seg_ele_pt])
3742  {
3743  last_seg_node_pt = last_seg_ele_pt->node_pt(0);
3744  }
3745 
3746  // Get the coordinates for the first and last segment's node
3747  for (unsigned i = 0; i < 2; i++)
3748  {
3749  first_seg_coord[i] = first_seg_node_pt->x(i);
3750  last_seg_coord[i] = last_seg_node_pt->x(i);
3751  }
3752 
3753  // -----------------------------------------------------------------
3754  // Copy the info. if the segment is inverted
3755  Boundary_segment_inverted[b].push_back(segment_inverted[is]);
3756 
3757  // Check if the segment is inverted, if that is the case then invert
3758  // the first and last seg. coordinates
3759  if (!segment_inverted[is])
3760  {
3761  // Store the initial and final coordinates that will help to
3762  // identify the segments in the new meshes created from this one
3763  Boundary_segment_initial_coordinate[b].push_back(first_seg_coord);
3764  Boundary_segment_final_coordinate[b].push_back(last_seg_coord);
3765  }
3766  else
3767  {
3768  // Store the initial and final coordinates that will help to
3769  // identify the segments in the new meshes created from this one
3770  // Invert the initial and final coordinates
3771  Boundary_segment_initial_coordinate[b].push_back(last_seg_coord);
3772  Boundary_segment_final_coordinate[b].push_back(first_seg_coord);
3773  }
3774 
3775  // Now assign initial and final zeta boundary coordinates for each
3776  // segment
3777  // -----------------------------------------------------------------
3778  // If there is a geom object then
3779  if (boundary_geom_object_pt(b) != 0)
3780  {
3781  // Store the initial and final zeta for the current segments (we
3782  // got this when we assigned arclength to the segments in the
3783  // current processor)
3784  if (segment_inverted[is])
3785  {
3786  Boundary_segment_initial_zeta[b].push_back(final_zeta_segment[is]);
3787  Boundary_segment_final_zeta[b].push_back(initial_zeta_segment[is]);
3788  }
3789  else
3790  {
3791  Boundary_segment_initial_zeta[b].push_back(initial_zeta_segment[is]);
3792  Boundary_segment_final_zeta[b].push_back(final_zeta_segment[is]);
3793  }
3794  } // if (boundary_geom_object_pt(b)!=0)
3795  else
3796  {
3797  // Store the initial arclength and vertices number for the
3798  // current segment
3799  Boundary_segment_initial_arclength[b].push_back(
3800  initial_segment_arclength[is]);
3801 
3802  Boundary_segment_final_arclength[b].push_back(
3803  initial_segment_arclength[is] + segment_arclength[is]);
3804 
3805  } // else if (boundary_geom_object_pt(b)!=0)
3806 
3807  } // // for (is < nsegments)
3808 
3809  // Get the number of segments from the sets of nodes
3810 #ifdef PARANOID
3811  if (segment_all_nodes_pt.size() != nsegments)
3812  {
3813  std::ostringstream error_message;
3814  error_message << "The number of segments (" << nsegments
3815  << ") and the number of "
3816  << "set of nodes (" << segment_all_nodes_pt.size()
3817  << ") representing\n"
3818  << "the\nsegments is different!!!\n\n";
3819  throw OomphLibError(error_message.str(),
3820  "TriangleMesh::compute_boundary_segments_"
3821  "connectivity_and_initial_zeta_values()",
3822  OOMPH_EXCEPTION_LOCATION);
3823  }
3824 #endif
3825 
3826  // The nodes have been assigned arc-length coordinates from one end
3827  // or the other of the connected segment.
3828 
3829  // -----------------------------------------------------------------
3830  // If mesh is distributed get the info. regarding the initial and
3831  // final nodes coordinates on the boundary, same as the zeta
3832  // boundary values for those nodes
3833 
3834  // Storage for the coordinates of the first and last nodes on the
3835  // boundary
3836  Vector<double> first_coordinate(2);
3837  Vector<double> last_coordinate(2);
3838 
3839  // Storage for the zeta coordinate of the first and last nodes on
3840  // the boundary
3841  Vector<double> first_node_zeta_coordinate(1, 0.0);
3842  Vector<double> last_node_zeta_coordinate(1, 0.0);
3843 
3844  // Send three data to all processors, the x[0], x[1] coordinate and
3845  // the zeta coordinate
3846  const unsigned ndtotal_data = 3;
3847  Vector<double> flat_packed_double_data_initial_seg(ndtotal_data);
3848 
3849  // If the mesh is distributed then check if this processor has the
3850  // initial segment
3851  if (my_rank == proc_with_initial_seg)
3852  {
3853  // Stores the firts element of the segment
3854  FiniteElement* first_ele_pt = 0;
3855  // Stores the first node of the boundary
3856  Node* first_node_pt = 0;
3857  // Check if the segment is inverted
3858  if (!segment_inverted[initial_segment])
3859  {
3860  // Get access to the first element on the segment marked as
3861  // initial
3862  first_ele_pt = segment_sorted_ele_pt[initial_segment].front();
3863 
3864  // Number of nodes
3865  const unsigned nnod = first_ele_pt->nnode();
3866 
3867  // Get the first node of the current segment
3868  first_node_pt = first_ele_pt->node_pt(0);
3869  if (is_inverted[first_ele_pt])
3870  {
3871  first_node_pt = first_ele_pt->node_pt(nnod - 1);
3872  }
3873  } // if (!segment_inverted[initial_segment])
3874  else
3875  {
3876  // Get access to the first element on the segment marked as
3877  // initial
3878  first_ele_pt = segment_sorted_ele_pt[initial_segment].back();
3879 
3880  // Number of nodes
3881  const unsigned nnod = first_ele_pt->nnode();
3882 
3883  // Get the first node of the current segment
3884  first_node_pt = first_ele_pt->node_pt(nnod - 1);
3885  if (is_inverted[first_ele_pt])
3886  {
3887  first_node_pt = first_ele_pt->node_pt(0);
3888  }
3889  } // else if (!segment_inverted[initial_segment])
3890 
3891  // Get the coordinates for the first node
3892  for (unsigned i = 0; i < 2; i++)
3893  {
3894  flat_packed_double_data_initial_seg[i] = first_node_pt->x(i);
3895  }
3896 
3897  // Get the zeta coordinates for the first node
3898  Vector<double> tmp_zeta(1);
3899  first_node_pt->get_coordinates_on_boundary(b, tmp_zeta);
3900 
3901  // If there is a geometric object associated to the boundary then
3902  // further process is necessary
3903  if (this->boundary_geom_object_pt(b) != 0)
3904  {
3905  // tmp_zeta[0] = this->boundary_coordinate_limits(b)[0];
3906  }
3907  else
3908  {
3909  // Check if the initial boundary coordinate is different from
3910  // zero, if that is the case then we need to set it to zero
3911  if (tmp_zeta[0] >= 1.0e-14)
3912  {
3913  tmp_zeta[0] = 0;
3914  }
3915  } // if (this->boundary_geom_object_pt(b)!=0)
3916 
3917  // Store the initial zeta value
3918  flat_packed_double_data_initial_seg[2] = tmp_zeta[0];
3919 
3920  } // if (my_rank == proc_with_initial_seg)
3921 
3922  // All processor receive the info. from the processor that has the
3923  // initial segment
3924  MPI_Bcast(&flat_packed_double_data_initial_seg[0],
3925  ndtotal_data,
3926  MPI_DOUBLE,
3927  proc_with_initial_seg,
3928  comm_pt->mpi_comm());
3929 
3930  // ... and all processor put that info. into the appropriate
3931  // storages
3932  for (unsigned i = 0; i < 2; i++)
3933  {
3934  first_coordinate[i] = flat_packed_double_data_initial_seg[i];
3935  }
3936  first_node_zeta_coordinate[0] = flat_packed_double_data_initial_seg[2];
3937 
3938  // -----------------------------------------------------------------
3939  // Send three data to all processors, the x[0], x[1] coordinate and
3940  // the zeta coordinate
3941  Vector<double> flat_packed_double_data_final_seg(ndtotal_data);
3942 
3943  // If the mesh is distributed then check if this processor has the
3944  // final segment
3945  if (my_rank == proc_with_final_seg)
3946  {
3947  // Get access to the last element on the segment
3948  FiniteElement* last_ele_pt = 0;
3949 
3950  // Get the last node of the current segment
3951  Node* last_node_pt = 0;
3952 
3953  // Check if the segment is inverted
3954  if (!segment_inverted[final_segment])
3955  {
3956  // Get access to the last element on the segment marked as
3957  // final
3958  last_ele_pt = segment_sorted_ele_pt[final_segment].back();
3959 
3960  // Number of nodes
3961  const unsigned nnod = last_ele_pt->nnode();
3962 
3963  // Get the last node of the current segment
3964  last_node_pt = last_ele_pt->node_pt(nnod - 1);
3965  if (is_inverted[last_ele_pt])
3966  {
3967  last_node_pt = last_ele_pt->node_pt(0);
3968  }
3969  } // if (!segment_inverted[final_segment])
3970  else
3971  {
3972  // Get access to the first element on the segment marked as
3973  // initial
3974  last_ele_pt = segment_sorted_ele_pt[final_segment].front();
3975 
3976  // Number of nodes
3977  const unsigned nnod = last_ele_pt->nnode();
3978 
3979  // Get the first node of the current segment
3980  last_node_pt = last_ele_pt->node_pt(0);
3981  if (is_inverted[last_ele_pt])
3982  {
3983  last_node_pt = last_ele_pt->node_pt(nnod - 1);
3984  }
3985  } // if (!segment_inverted[final_segment])
3986 
3987  // Get the coordinates for the last node
3988  for (unsigned i = 0; i < 2; i++)
3989  {
3990  flat_packed_double_data_final_seg[i] = last_node_pt->x(i);
3991  }
3992 
3993  // Get the zeta coordinates for the last node
3994  Vector<double> tmp_zeta(1);
3995  last_node_pt->get_coordinates_on_boundary(b, tmp_zeta);
3996 
3997  // If there is not a geometric object associated to the boundary
3998  // then further process is required
3999  if (this->boundary_geom_object_pt(b) != 0)
4000  {
4001  // Do nothing
4002  } // if (this->boundary_geom_object_pt(b)!=0)
4003  else
4004  {
4005  // Check if the final boundary coordinate is different from
4006  // the boundary arclength, if that is the case then we need
4007  // to set it to the accumulated arclength
4008  if (std::fabs(tmp_zeta[0] - root_accumulated_arclength) >= 1.0e-14)
4009  {
4010  tmp_zeta[0] = root_accumulated_arclength;
4011  }
4012  } // else if (this->boundary_geom_object_pt(b)!=0)
4013 
4014  // Store the final zeta value
4015  flat_packed_double_data_final_seg[2] = tmp_zeta[0];
4016 
4017  } // if (my_rank == proc_with_final_seg)
4018 
4019  // All processor receive the info. from the processor that has the
4020  // final segment
4021  MPI_Bcast(&flat_packed_double_data_final_seg[0],
4022  ndtotal_data,
4023  MPI_DOUBLE,
4024  proc_with_final_seg,
4025  comm_pt->mpi_comm());
4026 
4027  // All processor receive the info. from the processor that has the
4028  // final segment
4029  for (unsigned i = 0; i < 2; i++)
4030  {
4031  last_coordinate[i] = flat_packed_double_data_final_seg[i];
4032  }
4033  last_node_zeta_coordinate[0] = flat_packed_double_data_final_seg[2];
4034 
4035  // -----------------------------------------------------------------
4036  // Copy the values to the permanent storage
4037  Boundary_initial_coordinate[b] = first_coordinate;
4038  Boundary_final_coordinate[b] = last_coordinate;
4039 
4040  Boundary_initial_zeta_coordinate[b] = first_node_zeta_coordinate;
4041  Boundary_final_zeta_coordinate[b] = last_node_zeta_coordinate;
4042 
4043  // If we are dealing with an internal boundary then re-assign the
4044  // initial and final zeta values for the segments
4045  if (is_internal_boundary)
4046  {
4047  // Only re-assign zeta values if there are at least one nonhalo
4048  // segment, if all the possible segments are halo then the
4049  // synchronisation method will be in charge of assigning the
4050  // correct boundary coordinates
4051  if (nsegments > 0)
4052  {
4053  // Call the following method to re-construct the segments but
4054  // using only the nonhalo elements, therefore the boundary
4055  // coordinates need to be re-assigned
4056  re_assign_initial_zeta_values_for_internal_boundary(
4057  b, segment_sorted_nonhalo_ele_pt, is_inverted);
4058  }
4059 
4060  } // if (is_internal_boundary)
4061 
4062  // Now identify the boundary segments
4063  if (nsegments > 0)
4064  {
4065  // Identify the boundary segments in the current mesh
4066  // identify_boundary_segments_and_assign_initial_zeta_values(
4067  // b, all_face_ele_pt, is_internal_boundary, face_to_bulk_element_pt);
4068  identify_boundary_segments_and_assign_initial_zeta_values(b, this);
4069  } // if (nsegments > 0)
4070 
4071  // Clean all the created face elements
4072  for (unsigned i = 0; i < n_all_face_ele; i++)
4073  {
4074  delete all_face_ele_pt[i];
4075  all_face_ele_pt[i] = 0;
4076  }
4077  }
4078 
4079  //======================================================================
4080  /// Re-assign the boundary segments initial zeta (arclength)
4081  /// for those internal boundaries that were splited during the
4082  /// distribution process. Those boundaries that have one face element
4083  /// at each side of the boundary. Here we create the segments only
4084  /// with the nonhalo elements, therefore the boundary coordinates
4085  /// need to be re-assigned to be passed to the new meshes
4086  //======================================================================
4087  template<class ELEMENT>
4090  const unsigned& b,
4091  Vector<std::list<FiniteElement*>>& old_segment_sorted_ele_pt,
4092  std::map<FiniteElement*, bool>& old_is_inverted)
4093  {
4094  // ------------------------------------------------------------------
4095  // First: Get the face elements associated with the current boundary
4096  // Only include nonhalo face elements
4097  // ------------------------------------------------------------------
4098  // Temporary storage for face elements
4099  Vector<FiniteElement*> face_el_pt;
4100 
4101  // Temporary storage for the number of elements adjacent to the
4102  // boundary
4103  unsigned nele = 0;
4104 
4105  // Temporary storage for elements adjacent to the boundary that have
4106  // a common edge (related with internal boundaries)
4107  unsigned n_repeated_ele = 0;
4108 
4109  const unsigned n_regions = this->nregion();
4110 
4111  // Temporary storage for already done nodes
4112  Vector<std::pair<Node*, Node*>> done_nodes_pt;
4113 
4114  // If there is more than one region then only use boundary
4115  // coordinates from the bulk side (region 0)
4116  if (n_regions > 1)
4117  {
4118  for (unsigned rr = 0; rr < n_regions; rr++)
4119  {
4120  const unsigned region_id =
4121  static_cast<unsigned>(this->Region_attribute[rr]);
4122 
4123  // Loop over all elements on boundaries in region i_r
4124  const unsigned nel_in_region =
4125  this->nboundary_element_in_region(b, region_id);
4126 
4127  unsigned nel_repetead_in_region = 0;
4128 
4129  // Only bother to do anything else, if there are elements
4130  // associated with the boundary and the current region
4131  if (nel_in_region > 0)
4132  {
4133  bool repeated = false;
4134 
4135  // Loop over the bulk elements adjacent to boundary b
4136  for (unsigned e = 0; e < nel_in_region; e++)
4137  {
4138  // Get pointer to the bulk element that is adjacent to
4139  // boundary b
4140  FiniteElement* bulk_elem_pt =
4141  this->boundary_element_in_region_pt(b, region_id, e);
4142 
4143  // Remember only work with non halo elements
4144  if (bulk_elem_pt->is_halo())
4145  {
4146  n_repeated_ele++;
4147  continue;
4148  }
4149 
4150  // Find the index of the face of element e along boundary b
4151  int face_index =
4152  this->face_index_at_boundary_in_region(b, region_id, e);
4153 
4154  // Before adding the new element we need to be sure that the
4155  // edge that this element represent has not been already
4156  // added
4157  FiniteElement* tmp_ele_pt =
4158  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
4159 
4160  const unsigned n_nodes = tmp_ele_pt->nnode();
4161 
4162  std::pair<Node*, Node*> tmp_pair = std::make_pair(
4163  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
4164 
4165  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
4166  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
4167 
4168  // Search for repeated nodes
4169  const unsigned repeated_nodes_size = done_nodes_pt.size();
4170  for (unsigned l = 0; l < repeated_nodes_size; l++)
4171  {
4172  if (tmp_pair == done_nodes_pt[l] ||
4173  tmp_pair_inverse == done_nodes_pt[l])
4174  {
4175  nel_repetead_in_region++;
4176  repeated = true;
4177  break;
4178  }
4179  }
4180 
4181  // Create new face element
4182  if (!repeated)
4183  {
4184  // Add the pair of nodes (edge) to the node dones
4185  done_nodes_pt.push_back(tmp_pair);
4186  // Add the element to the face elements
4187  face_el_pt.push_back(tmp_ele_pt);
4188  }
4189  else
4190  {
4191  // Clean up
4192  delete tmp_ele_pt;
4193  tmp_ele_pt = 0;
4194  }
4195 
4196  // Re-start
4197  repeated = false;
4198 
4199  } // for nel
4200 
4201  nele += nel_in_region;
4202 
4203  n_repeated_ele += nel_repetead_in_region;
4204 
4205  } // if (nel_in_region > 0)
4206  } // for (rr < n_regions)
4207  } // if (n_regions > 1)
4208  // Otherwise it's just the normal boundary functions
4209  else
4210  {
4211  // Loop over all elements on boundaries
4212  nele = this->nboundary_element(b);
4213 
4214  // Only bother to do anything else, if there are elements
4215  if (nele > 0)
4216  {
4217  // Check for repeated ones
4218  bool repeated = false;
4219 
4220  // Loop over the bulk elements adjacent to boundary b
4221  for (unsigned e = 0; e < nele; e++)
4222  {
4223  // Get pointer to the bulk element that is adjacent to
4224  // boundary b
4225  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
4226 
4227  // Skip the halo elements, they are not included
4228  if (bulk_elem_pt->is_halo())
4229  {
4230  n_repeated_ele++;
4231  continue;
4232  }
4233 
4234  // Find the index of the face of element e along boundary b
4235  int face_index = this->face_index_at_boundary(b, e);
4236 
4237  // Before adding the new element we need to be sure that the
4238  // edge that this element represents has not been already
4239  // added (only applies for internal boundaries)
4240  FiniteElement* tmp_ele_pt =
4241  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
4242 
4243  const unsigned n_nodes = tmp_ele_pt->nnode();
4244 
4245  std::pair<Node*, Node*> tmp_pair = std::make_pair(
4246  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
4247 
4248  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
4249  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
4250 
4251  // Search for repeated nodes
4252  const unsigned repeated_nodes_size = done_nodes_pt.size();
4253  for (unsigned l = 0; l < repeated_nodes_size; l++)
4254  {
4255  if (tmp_pair == done_nodes_pt[l] ||
4256  tmp_pair_inverse == done_nodes_pt[l])
4257  {
4258  // Increase the number of repeated elements
4259  n_repeated_ele++;
4260  // Mark the element as repeated
4261  repeated = true;
4262  break;
4263  }
4264  }
4265 
4266  // Create new face element
4267  if (!repeated)
4268  {
4269  // Add the pair of nodes (edge) to the node dones
4270  done_nodes_pt.push_back(tmp_pair);
4271  // Add the element to the face elements
4272  face_el_pt.push_back(tmp_ele_pt);
4273  }
4274  else
4275  {
4276  // Free the repeated bulk element!!
4277  delete tmp_ele_pt;
4278  tmp_ele_pt = 0;
4279  }
4280 
4281  // Re-start
4282  repeated = false;
4283 
4284  } // for (e < nel)
4285  } // if (nel > 0)
4286 
4287  } // else (n_regions > 1)
4288 
4289  // Do not consider the repeated elements
4290  nele -= n_repeated_ele;
4291 
4292 #ifdef PARANOID
4293  if (nele != face_el_pt.size())
4294  {
4295  std::ostringstream error_message;
4296  error_message
4297  << "The independet counting of face elements (" << nele << ") for "
4298  << "boundary (" << b << ") is different\n"
4299  << "from the real number of face elements in the container ("
4300  << face_el_pt.size() << ")\n";
4301  //<< "Possible memory leak\n"
4302  throw OomphLibError(
4303  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4304  }
4305 #endif
4306 
4307  // ----------------------------------------------------------------
4308  // Second: Sort the face elements, only consider nonhalo elements
4309  // ----------------------------------------------------------------
4310 
4311  // Get the total number of nonhalo face elements
4312  const unsigned nnon_halo_face_elements = face_el_pt.size();
4313 
4314  // The vector of list to store the "segments" that compound the
4315  // boundary (segments may appear only in a distributed mesh)
4316  Vector<std::list<FiniteElement*>> segment_sorted_ele_pt;
4317 
4318  // Number of already sorted face elements
4319  unsigned nsorted_face_elements = 0;
4320 
4321  // Keep track of who's done
4322  std::map<FiniteElement*, bool> done_el;
4323 
4324  // Keep track of which element is inverted
4325  std::map<FiniteElement*, bool> is_inverted;
4326 
4327  // Iterate until all possible segments have been created
4328  while (nsorted_face_elements < nnon_halo_face_elements)
4329  {
4330  // The ordered list of face elements (in a distributed mesh a
4331  // collection of contiguous face elements define a segment)
4332  std::list<FiniteElement*> sorted_el_pt;
4333 
4334 #ifdef PARANOID
4335  // Select an initial element for the segment
4336  bool found_initial_face_element = false;
4337 #endif
4338 
4339  FiniteElement* ele_face_pt = 0;
4340 
4341  unsigned iface = 0;
4342  for (iface = 0; iface < nele; iface++)
4343  {
4344  ele_face_pt = face_el_pt[iface];
4345  // If not done then take it as initial face element
4346  if (!done_el[ele_face_pt])
4347  {
4348 #ifdef PARANOID
4349  // Mark as found the root face element
4350  found_initial_face_element = true;
4351 #endif
4352  // Increase the number of sorted face elements
4353  nsorted_face_elements++;
4354  // Increase the counter to mark the position of the next
4355  // element number
4356  iface++;
4357  // Add the face element in the list of sorted face elements
4358  sorted_el_pt.push_back(ele_face_pt);
4359  // Mark as done
4360  done_el[ele_face_pt] = true;
4361  break;
4362  } // if (!done_el[ele_face_pt])
4363  } // for (iface < nele)
4364 
4365 #ifdef PARANOID
4366  if (!found_initial_face_element)
4367  {
4368  std::ostringstream error_message;
4369  error_message
4370  << "Could not find an initial face element for the current segment\n";
4371  throw OomphLibError(
4372  error_message.str(),
4373  "TriangleMesh::re_assign_initial_zeta_values_for_internal_boundary()",
4374  OOMPH_EXCEPTION_LOCATION);
4375  }
4376 #endif
4377 
4378  // Number of nodes
4379  const unsigned nnod = ele_face_pt->nnode();
4380 
4381  // Left and rightmost nodes (the left and right nodes of the
4382  // current face element)
4383  Node* left_node_pt = ele_face_pt->node_pt(0);
4384  Node* right_node_pt = ele_face_pt->node_pt(nnod - 1);
4385 
4386  // Continue iterating if a new face element has been added to the
4387  // list
4388  bool face_element_added = false;
4389 
4390  // While a new face element has been added to the set of sorted
4391  // face elements then re-iterate
4392  do
4393  {
4394  // Start from the next face element since we have already added
4395  // the previous one as the initial face element (any previous
4396  // face element had to be added on previous iterations)
4397  for (unsigned iiface = iface; iiface < nele; iiface++)
4398  {
4399  // Re-start flag
4400  face_element_added = false;
4401 
4402  // Get the candidate element
4403  ele_face_pt = face_el_pt[iiface];
4404 
4405  // Check that the candidate element has not been done and is
4406  // not a halo element
4407  if (!(done_el[ele_face_pt]))
4408  {
4409  // Get the left and right nodes of the current element
4410  Node* local_left_node_pt = ele_face_pt->node_pt(0);
4411  Node* local_right_node_pt = ele_face_pt->node_pt(nnod - 1);
4412 
4413  // New element fits at the left of segment and is not inverted
4414  if (left_node_pt == local_right_node_pt)
4415  {
4416  left_node_pt = local_left_node_pt;
4417  sorted_el_pt.push_front(ele_face_pt);
4418  is_inverted[ele_face_pt] = false;
4419  face_element_added = true;
4420  }
4421  // New element fits at the left of segment and is inverted
4422  else if (left_node_pt == local_left_node_pt)
4423  {
4424  left_node_pt = local_right_node_pt;
4425  sorted_el_pt.push_front(ele_face_pt);
4426  is_inverted[ele_face_pt] = true;
4427  face_element_added = true;
4428  }
4429  // New element fits on the right of segment and is not inverted
4430  else if (right_node_pt == local_left_node_pt)
4431  {
4432  right_node_pt = local_right_node_pt;
4433  sorted_el_pt.push_back(ele_face_pt);
4434  is_inverted[ele_face_pt] = false;
4435  face_element_added = true;
4436  }
4437  // New element fits on the right of segment and is inverted
4438  else if (right_node_pt == local_right_node_pt)
4439  {
4440  right_node_pt = local_left_node_pt;
4441  sorted_el_pt.push_back(ele_face_pt);
4442  is_inverted[ele_face_pt] = true;
4443  face_element_added = true;
4444  }
4445 
4446  if (face_element_added)
4447  {
4448  done_el[ele_face_pt] = true;
4449  nsorted_face_elements++;
4450  break;
4451  } // if (face_element_added)
4452 
4453  } // if (!(done_el[ele_face_pt]))
4454 
4455  } // for (iiface<nnon_halo_face_element)
4456 
4457  } while (face_element_added &&
4458  (nsorted_face_elements < nnon_halo_face_elements));
4459 
4460  // Store the created segment in the vector of segments
4461  segment_sorted_ele_pt.push_back(sorted_el_pt);
4462 
4463  } // while(nsorted_face_elements < nnon_halo_face_elements);
4464 
4465  // --------------------------------------------------------------
4466  // Third: We have the face elements sorted, now assign boundary
4467  // coordinates to the nodes in the segments and compute the
4468  // arclength of the segment.
4469  // --------------------------------------------------------------
4470 
4471  // The number of segments in this processor
4472  const unsigned nsegments = segment_sorted_ele_pt.size();
4473 
4474 #ifdef PARANOID
4475  if (nnon_halo_face_elements > 0 && nsegments == 0)
4476  {
4477  std::ostringstream error_message;
4478  error_message
4479  << "The number of segments is zero, but the number of nonhalo\n"
4480  << "elements is: (" << nnon_halo_face_elements << ")\n";
4481  throw OomphLibError(
4482  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4483  } // if (nnon_halo_face_elements > 0 && nsegments == 0)
4484 #endif
4485 
4486  // Vector of sets that stores the nodes of each segment based on a
4487  // lexicographically order starting from the bottom left node of
4488  // each segment
4489  Vector<std::set<Node*>> segment_all_nodes_pt(nsegments);
4490 
4491  // Stores the nodes on each segment in the order they appear in the
4492  // face elements
4493  Vector<Vector<Node*>> sorted_segment_all_nodes_pt(nsegments);
4494 
4495  // Associate and arclength to each node on each segment of the
4496  // boundary, the nodes and therefore the arclength come in the same
4497  // order as the face elements
4498  Vector<Vector<double>> sorted_segment_node_arclength(nsegments);
4499 
4500  // The arclength of each segment in the current processor
4501  Vector<double> segment_arclength(nsegments);
4502 
4503  // The number of vertices of each segment
4504  Vector<unsigned> nvertices_per_segment(nsegments);
4505 
4506  // The initial zeta for the segment
4507  Vector<double> initial_zeta_segment(nsegments);
4508 
4509  // The final zeta for the segment
4510  Vector<double> final_zeta_segment(nsegments);
4511 
4512  // Go through all the segments and compute the LOCAL boundary
4513  // coordinates
4514  for (unsigned is = 0; is < nsegments; is++)
4515  {
4516 #ifdef PARANOID
4517  if (segment_sorted_ele_pt[is].size() == 0)
4518  {
4519  std::ostringstream error_message;
4520  error_message << "The (" << is << ")-th segment has no elements\n";
4521  throw OomphLibError(
4522  error_message.str(),
4523  "TriangleMesh::re_assign_initial_zeta_values_for_internal_boundary()",
4524  OOMPH_EXCEPTION_LOCATION);
4525  } // if (segment_sorted_ele_pt[is].size() == 0)
4526 #endif
4527 
4528  // Get access to the first element on the segment
4529  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
4530 
4531  // Number of nodes
4532  const unsigned nnod = first_ele_pt->nnode();
4533 
4534  // Get the first node of the current segment
4535  Node* first_node_pt = first_ele_pt->node_pt(0);
4536  if (is_inverted[first_ele_pt])
4537  {
4538  first_node_pt = first_ele_pt->node_pt(nnod - 1);
4539  }
4540 
4541  // Coordinates of left node
4542  double x_left = first_node_pt->x(0);
4543  double y_left = first_node_pt->x(1);
4544 
4545  // Initialise boundary coordinate (local boundary coordinate for
4546  // boundaries with more than one segment)
4547  Vector<double> zeta(1, 0.0);
4548 
4549  // If we have associated a GeomObject then it is not necessary
4550  // to compute the arclength, only read the values from the nodes at
4551  // the edges
4552  if (this->boundary_geom_object_pt(b) != 0)
4553  {
4554  first_node_pt->get_coordinates_on_boundary(b, zeta);
4555  initial_zeta_segment[is] = zeta[0];
4556 
4557  // Get access to the last element on the segment
4558  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
4559 
4560  // Get the last node of the current segment
4561  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
4562  if (is_inverted[last_ele_pt])
4563  {
4564  last_node_pt = last_ele_pt->node_pt(0);
4565  }
4566 
4567  last_node_pt->get_coordinates_on_boundary(b, zeta);
4568  final_zeta_segment[is] = zeta[0];
4569  }
4570 
4571  // Sort the nodes in the segment (lexicographically bottom left
4572  // node)
4573  std::set<Node*> local_nodes_pt;
4574  local_nodes_pt.insert(first_node_pt);
4575 
4576  // Associate and arclength to the sorted nodes
4577  Vector<double> sorted_node_arclength;
4578  sorted_node_arclength.push_back(0.0);
4579 
4580  // Sorts the nodes in the segments according their sorting in the
4581  // face elements
4582  Vector<Node*> sorted_nodes_pt;
4583  sorted_nodes_pt.push_back(first_node_pt);
4584 
4585  // Now loop over nodes in order
4586  for (std::list<FiniteElement*>::iterator it =
4587  segment_sorted_ele_pt[is].begin();
4588  it != segment_sorted_ele_pt[is].end();
4589  it++)
4590  {
4591  // Get the face element
4592  FiniteElement* el_pt = *it;
4593 
4594  // Start node and increment
4595  unsigned k_nod = 1;
4596  int nod_diff = 1;
4597  if (is_inverted[el_pt])
4598  {
4599  k_nod = nnod - 2;
4600  nod_diff = -1;
4601  }
4602 
4603  // Loop over nodes
4604  for (unsigned j = 1; j < nnod; j++)
4605  {
4606  Node* nod_pt = el_pt->node_pt(k_nod);
4607  k_nod += nod_diff;
4608 
4609  // Coordinates of right node
4610  double x_right = nod_pt->x(0);
4611  double y_right = nod_pt->x(1);
4612 
4613  // Increment boundary coordinate
4614  zeta[0] += sqrt((x_right - x_left) * (x_right - x_left) +
4615  (y_right - y_left) * (y_right - y_left));
4616 
4617  // When we have a GeomObject associated to the boundary we already
4618  // know the zeta values for the nodes, there is no need to compute
4619  // the arclength
4620  if (this->boundary_geom_object_pt(b) == 0)
4621  {
4622  // Set boundary coordinate
4623  // nod_pt->set_coordinates_on_boundary(b, zeta);
4624  }
4625 
4626  // Increment reference coordinate
4627  x_left = x_right;
4628  y_left = y_right;
4629 
4630  // Get lexicographically bottom left node but only
4631  // use vertex nodes as candidates
4632  local_nodes_pt.insert(nod_pt);
4633 
4634  // Associate the arclength for the current node
4635  sorted_node_arclength.push_back(zeta[0]);
4636 
4637  // Store the node in the sorted nodes storage
4638  sorted_nodes_pt.push_back(nod_pt);
4639 
4640  } // for (j < nnod)
4641 
4642  } // iterator over the elements in the segment
4643 
4644  // Info. to be passed to the other processors
4645  // The initial arclength for the segment that goes after this depends
4646  // on the current segment arclength
4647  segment_arclength[is] = zeta[0];
4648 
4649  // Info. to be passed to the other processors
4650  // The initial vertex number for the segment that goes after this
4651  // depends on the current sement vertices number
4652  nvertices_per_segment[is] = local_nodes_pt.size();
4653 
4654  // Add the nodes for the corresponding segment in the container
4655  segment_all_nodes_pt[is] = local_nodes_pt;
4656 
4657  // Add the arclengths to the nodes in the segment
4658  sorted_segment_node_arclength[is] = sorted_node_arclength;
4659 
4660  // Add the sorted nodes to the storage
4661  sorted_segment_all_nodes_pt[is] = sorted_nodes_pt;
4662 
4663  // The attaching of the halo elements at both sides of the segments is
4664  // performed only if segments connectivity needs to be computed
4665 
4666  } // for (is < nsegments)
4667 
4668  // ------------------------------------------------------------------
4669  // Fourth: Now we have the segments sorted, with arclength and with
4670  // LOCAL boundary coordinates assigned to the nodes. Identify the
4671  // nodes on the segments with the input segments and re-assign all
4672  // the info. related with the identification of segments
4673  // ------------------------------------------------------------------
4674 
4675  // Get the number of segments for the old sorted segments
4676  const unsigned old_nsegments = old_segment_sorted_ele_pt.size();
4677 
4678  // ------------------------------------------------------------------
4679  // Copy the old info. in temporary storages
4680  Vector<unsigned> old_boundary_segment_inverted(old_nsegments);
4681 
4682  Vector<Vector<double>> old_boundary_segment_initial_coordinate(
4683  old_nsegments);
4684  Vector<Vector<double>> old_boundary_segment_final_coordinate(old_nsegments);
4685 
4686  Vector<double> old_boundary_segment_initial_zeta(old_nsegments);
4687  Vector<double> old_boundary_segment_final_zeta(old_nsegments);
4688 
4689  Vector<double> old_boundary_segment_initial_arclength(old_nsegments);
4690  Vector<double> old_boundary_segment_final_arclength(old_nsegments);
4691 
4692  // Back-up the information
4693  for (unsigned old_is = 0; old_is < old_nsegments; old_is++)
4694  {
4695  old_boundary_segment_inverted[old_is] =
4696  boundary_segment_inverted(b)[old_is];
4697 
4698  old_boundary_segment_initial_coordinate[old_is].resize(2);
4699  old_boundary_segment_final_coordinate[old_is].resize(2);
4700  for (unsigned i = 0; i < 2; i++)
4701  {
4702  old_boundary_segment_initial_coordinate[old_is][i] =
4703  boundary_segment_initial_coordinate(b)[old_is][i];
4704 
4705  old_boundary_segment_final_coordinate[old_is][i] =
4706  boundary_segment_final_coordinate(b)[old_is][i];
4707  }
4708 
4709  // Check if the boundary has an associated GeomObject
4710  if (this->boundary_geom_object_pt(b) != 0)
4711  {
4712  old_boundary_segment_initial_zeta[old_is] =
4713  boundary_segment_initial_zeta(b)[old_is];
4714 
4715  old_boundary_segment_final_zeta[old_is] =
4716  boundary_segment_final_zeta(b)[old_is];
4717 
4718  } // if (this->boundary_geom_object_pt(b)!=0)
4719  else
4720  {
4721  old_boundary_segment_initial_arclength[old_is] =
4722  boundary_segment_initial_arclength(b)[old_is];
4723 
4724  old_boundary_segment_final_arclength[old_is] =
4725  boundary_segment_final_arclength(b)[old_is];
4726 
4727  } // else if (this->boundary_geom_object_pt(b)!=0)
4728 
4729  } // for (old_is < old_nsegments)
4730 
4731  // ------------------------------------------------------------------
4732  // Now clear the original storages
4733  Boundary_segment_inverted[b].clear();
4734  Boundary_segment_initial_coordinate[b].clear();
4735  Boundary_segment_final_coordinate[b].clear();
4736 
4737  Boundary_segment_initial_zeta[b].clear();
4738  Boundary_segment_final_zeta[b].clear();
4739 
4740  Boundary_segment_initial_arclength[b].clear();
4741  Boundary_segment_final_arclength[b].clear();
4742  // ------------------------------------------------------------------
4743  // .. and resize the storages for the new number of segments
4744  Boundary_segment_inverted[b].resize(nsegments);
4745  Boundary_segment_initial_coordinate[b].resize(nsegments);
4746  Boundary_segment_final_coordinate[b].resize(nsegments);
4747 
4748  // Check if the boundary has an associated GeomObject
4749  if (this->boundary_geom_object_pt(b) != 0)
4750  {
4751  Boundary_segment_initial_zeta[b].resize(nsegments);
4752  Boundary_segment_final_zeta[b].resize(nsegments);
4753  }
4754  else
4755  {
4756  Boundary_segment_initial_arclength[b].resize(nsegments);
4757  Boundary_segment_final_arclength[b].resize(nsegments);
4758  }
4759  // ------------------------------------------------------------------
4760  // map to know if the new segment has been re-assigned the info.
4761  std::map<unsigned, bool> done_segment;
4762 
4763  // Count the number of re-assigned segments with the new values
4764  unsigned re_assigned_segments = 0;
4765 
4766  // Go through all the old segments (the input segments)
4767  for (unsigned old_is = 0; old_is < old_nsegments; old_is++)
4768  {
4769  // Get the first and last zeta values for the current segment
4770  const double old_initial_arclength =
4771  old_boundary_segment_initial_arclength[old_is];
4772  const double old_final_arclength =
4773  old_boundary_segment_final_arclength[old_is];
4774  // Get the "is inverted" segment information
4775  const unsigned old_inverted_segment =
4776  old_boundary_segment_inverted[old_is];
4777 
4778  // Check if the boundary coordinates in the segment go in
4779  // increasing or decreasing order
4780  bool old_increasing_order = false;
4781  if (old_initial_arclength < old_final_arclength)
4782  {
4783  old_increasing_order = true;
4784  }
4785 
4786  // Now get the first and last node of the current segment
4787  // Get the first element
4788  FiniteElement* first_old_seg_ele_pt =
4789  old_segment_sorted_ele_pt[old_is].front();
4790 
4791  // Number of nodes
4792  const unsigned nnod = first_old_seg_ele_pt->nnode();
4793 
4794  // Get the first node of the current segment
4795  Node* first_old_seg_node_pt = first_old_seg_ele_pt->node_pt(0);
4796  if (old_is_inverted[first_old_seg_ele_pt])
4797  {
4798  first_old_seg_node_pt = first_old_seg_ele_pt->node_pt(nnod - 1);
4799  }
4800 
4801  // Get access to the last element on the segment
4802  FiniteElement* last_old_seg_ele_pt =
4803  old_segment_sorted_ele_pt[old_is].back();
4804 
4805  // Get the last node of the current segment
4806  Node* last_old_seg_node_pt = last_old_seg_ele_pt->node_pt(nnod - 1);
4807  if (old_is_inverted[last_old_seg_ele_pt])
4808  {
4809  last_old_seg_node_pt = last_old_seg_ele_pt->node_pt(0);
4810  }
4811  // Check if the segment is inverted, if that is the case then
4812  // also invert the nodes
4813  if (old_inverted_segment)
4814  {
4815  Node* temp_node_pt = first_old_seg_node_pt;
4816  first_old_seg_node_pt = last_old_seg_node_pt;
4817  last_old_seg_node_pt = temp_node_pt;
4818  }
4819 
4820  // We have the first and last node of the old segment (input
4821  // segment), now identify in which segment, of those with only
4822  // nonhalo face elements, they are
4823  for (unsigned is = 0; is < nsegments; is++)
4824  {
4825  if (!done_segment[is])
4826  {
4827  // Go through the nodes of the current segment and try to find
4828  // the old nodes
4829  bool found_first_old_seg_node = false;
4830  bool found_last_old_seg_node = false;
4831  bool same_order = false;
4832 
4833  // Get the first node of the current segment
4834  FiniteElement* first_seg_ele_pt = segment_sorted_ele_pt[is].front();
4835  Node* first_seg_node_pt = first_seg_ele_pt->node_pt(0);
4836  if (is_inverted[first_seg_ele_pt])
4837  {
4838  first_seg_node_pt = first_seg_ele_pt->node_pt(nnod - 1);
4839  }
4840 
4841  // Get the arclength for the first node
4842  const double segment_first_node_zeta =
4843  sorted_segment_node_arclength[is][0];
4844 
4845  // Get the node coordinates for the first node
4846  Vector<double> first_node_coord(2);
4847  for (unsigned i = 0; i < 2; i++)
4848  {
4849  first_node_coord[i] = first_seg_node_pt->x(i);
4850  }
4851 
4852  // Get the last node of the current segment
4853  FiniteElement* last_seg_ele_pt = segment_sorted_ele_pt[is].back();
4854  Node* last_seg_node_pt = last_seg_ele_pt->node_pt(nnod - 1);
4855  if (is_inverted[last_seg_ele_pt])
4856  {
4857  last_seg_node_pt = last_seg_ele_pt->node_pt(0);
4858  }
4859 
4860  // Get the arclength for the last node
4861  const double segment_final_node_zeta = segment_arclength[is];
4862 
4863  // Get the node coordinates for the last node
4864  Vector<double> last_node_coord(2);
4865  for (unsigned i = 0; i < 2; i++)
4866  {
4867  last_node_coord[i] = last_seg_node_pt->x(i);
4868  }
4869 
4870  // Temporary storage for the nodes of the current segment
4871  Vector<Node*> segment_node_pt = sorted_segment_all_nodes_pt[is];
4872  // Get the number of nodes in the segment
4873  const unsigned nsegment_node = segment_node_pt.size();
4874  for (unsigned in = 0; in < nsegment_node; in++)
4875  {
4876  Node* current_node_pt = segment_node_pt[in];
4877  if (!found_first_old_seg_node &&
4878  first_old_seg_node_pt == current_node_pt)
4879  {
4880  // Get the arclength assigned to the node on the old
4881  // segment
4882  const double current_node_zeta =
4883  sorted_segment_node_arclength[is][in];
4884 
4885  // Now check if the new segment has the same orientation
4886  // as the old one
4887  if (!found_last_old_seg_node) // has the same orientation
4888  {
4889  // Re-assign the first node coordinates
4890  Boundary_segment_initial_coordinate[b][is] = first_node_coord;
4891 
4892  // Check if the boundary has an associated GeomObject
4893  if (this->boundary_geom_object_pt(b) != 0)
4894  {
4895  // Assign the zeta values if the current segment has the
4896  // nodes of the old one
4897 
4898  // If we are in the same order then pass the values as
4899  // they are
4900  Boundary_segment_initial_zeta[b][is] =
4901  initial_zeta_segment[is];
4902 
4903  } // if (this->boundary_geom_object_pt(b)!=0)
4904  else
4905  {
4906  // Get the distance to the first node
4907  const double distance =
4908  std::fabs(current_node_zeta - segment_first_node_zeta);
4909 
4910  double new_initial_arclength = old_initial_arclength;
4911 
4912  // Now check if the zeta values are in increasing order
4913  if (old_increasing_order)
4914  {
4915  // Substract the distance
4916  new_initial_arclength -= distance;
4917  }
4918  else
4919  {
4920  // Add the distance
4921  new_initial_arclength += distance;
4922  }
4923 
4924  // Re-assign the initial arclength for the current segment
4925  Boundary_segment_initial_arclength[b][is] =
4926  new_initial_arclength;
4927 
4928  } // else if (this->boundary_geom_object_pt(b)!=0)
4929  } // if (!found_last_old_seg_node)
4930  else // has different orientation
4931  {
4932  // Re-assign the first node coordinates
4933  Boundary_segment_initial_coordinate[b][is] = last_node_coord;
4934 
4935  // Check if the boundary has an associated GeomObject
4936  if (this->boundary_geom_object_pt(b) != 0)
4937  {
4938  // Assign the zeta values if the current segment has the
4939  // nodes of the old one
4940 
4941  // Not the same order, we need to copy the zeta values
4942  // from the other end, the inverted flag is changed at
4943  // the end. Copy the value from the final end
4944  Boundary_segment_initial_zeta[b][is] = final_zeta_segment[is];
4945 
4946  } // if (this->boundary_geom_object_pt(b)!=0)
4947  else
4948  {
4949  // Get the distance to the final node
4950  const double distance =
4951  std::fabs(current_node_zeta - segment_final_node_zeta);
4952 
4953  double new_initial_arclength = old_initial_arclength;
4954 
4955  // Now check if the zeta values are in increasing order
4956  if (old_increasing_order)
4957  {
4958  // Substract the distance
4959  new_initial_arclength -= distance;
4960  }
4961  else
4962  {
4963  // Add the distance
4964  new_initial_arclength += distance;
4965  }
4966 
4967  // Re-assign the initial arclength for the current segment
4968  Boundary_segment_initial_arclength[b][is] =
4969  new_initial_arclength;
4970 
4971  } // else if (this->boundary_geom_object_pt(b)!=0)
4972  } // else if (!found_last_old_seg_node)
4973 
4974  // Mark as found the first node
4975  found_first_old_seg_node = true;
4976  }
4977  // if (!found_first_old_seg_node &&
4978  // first_old_seg_node_pt == current_node_pt)
4979 
4980  // If we found first the first node then the segments have
4981  // the same order
4982  if (found_first_old_seg_node && !found_last_old_seg_node)
4983  {
4984  same_order = true;
4985  }
4986 
4987  if (!found_last_old_seg_node &&
4988  last_old_seg_node_pt == current_node_pt)
4989  {
4990  // Get the boundary coordinates assigned to the node on
4991  // the old segment
4992  const double current_node_zeta =
4993  sorted_segment_node_arclength[is][in];
4994 
4995  // Now check if the new segment has the same orientation
4996  // as the old one
4997  if (found_first_old_seg_node) // has the same orientation
4998  {
4999  // Re-assign the last node coordinates
5000  Boundary_segment_final_coordinate[b][is] = last_node_coord;
5001 
5002  // Check if the boundary has an associated GeomObject
5003  if (this->boundary_geom_object_pt(b) != 0)
5004  {
5005  // Assign the zeta values if the current segment has the
5006  // nodes of the old one
5007 
5008  // If we are in the same order then pass the values as
5009  // they are
5010  Boundary_segment_final_zeta[b][is] = final_zeta_segment[is];
5011 
5012  } // if (this->boundary_geom_object_pt(b)!=0)
5013  else
5014  {
5015  // Get the distance to the last node
5016  const double distance =
5017  std::fabs(current_node_zeta - segment_final_node_zeta);
5018 
5019  double new_final_arclength = old_final_arclength;
5020 
5021  // Now check if the zeta values are in increasing order
5022  if (old_increasing_order)
5023  {
5024  // Add the distance
5025  new_final_arclength += distance;
5026  }
5027  else
5028  {
5029  // Substract the distance
5030  new_final_arclength -= distance;
5031  }
5032 
5033  // Re-assign the final arclength for the current segment
5034  Boundary_segment_final_arclength[b][is] = new_final_arclength;
5035 
5036  } // else if (this->boundary_geom_object_pt(b)!=0)
5037  } // if (found_first_old_seg_node)
5038  else
5039  {
5040  // Re-assign the last node coordinates
5041  Boundary_segment_final_coordinate[b][is] = first_node_coord;
5042 
5043  // Check if the boundary has an associated GeomObject
5044  if (this->boundary_geom_object_pt(b) != 0)
5045  {
5046  // Assign the zeta values if the current segment has the
5047  // nodes of the old one
5048 
5049  // Not the same order, we need to copy the zeta values
5050  // from the other end, the inverted flag is changed at
5051  // the end. Copy the value from the initial end
5052  Boundary_segment_final_zeta[b][is] = initial_zeta_segment[is];
5053 
5054  } // if (this->boundary_geom_object_pt(b)!=0)
5055  else
5056  {
5057  // Get the distance to the last node
5058  const double distance =
5059  std::fabs(current_node_zeta - segment_first_node_zeta);
5060 
5061  double new_final_arclength = old_final_arclength;
5062 
5063  // Now check if the zeta values are in increasing order
5064  if (old_increasing_order)
5065  {
5066  // Add the distance
5067  new_final_arclength += distance;
5068  }
5069  else
5070  {
5071  // Substract the distance
5072  new_final_arclength -= distance;
5073  }
5074 
5075  // Re-assign the final arclength for the current segment
5076  Boundary_segment_final_arclength[b][is] = new_final_arclength;
5077 
5078  } // else if (this->boundary_geom_object_pt(b)!=0)
5079  } // if (found_first_old_seg_node)
5080 
5081  // Mark as found the last node
5082  found_last_old_seg_node = true;
5083 
5084  } // if (!found_last_old_seg_node &&
5085  // last_old_seg_node_pt == current_node_pt)
5086 
5087  // If we found the last node first then the segments have
5088  // not the same order
5089  if (!found_first_old_seg_node && found_last_old_seg_node)
5090  {
5091  same_order = false;
5092  }
5093 
5094  if (found_first_old_seg_node && found_last_old_seg_node)
5095  {
5096  // Check if necessary to change the information that
5097  // states if a segment is inverted or not
5098  if (same_order)
5099  {
5100  Boundary_segment_inverted[b][is] = old_inverted_segment;
5101  }
5102  else
5103  {
5104  Boundary_segment_inverted[b][is] = !old_inverted_segment;
5105  }
5106 
5107  // Mark the segment as done
5108  done_segment[is] = true;
5109 
5110  // Increase the number of re-assigned segments
5111  re_assigned_segments++;
5112 
5113  // Break the for that look for the nodes in the segments
5114  break;
5115  }
5116 
5117  } // for (in < nsegment_node)
5118 
5119 #ifdef PARANOID
5120  if ((found_first_old_seg_node && !found_last_old_seg_node) ||
5121  (!found_first_old_seg_node && found_last_old_seg_node))
5122  {
5123  std::stringstream error_message;
5124  error_message
5125  << "Working with boundary (" << b << ").\nOnly the first node or "
5126  << "the last node of the old segment (" << old_is << ") was\n"
5127  << "found. Both, first and last node should have been found in "
5128  << "the same segment!!!.\n"
5129  << "Found first seg node:" << found_first_old_seg_node << "\n"
5130  << "Found last seg node:" << found_last_old_seg_node << "\n\n";
5131  throw OomphLibError(error_message.str(),
5132  "TriangleMesh::re_assign_initial_zeta_values_"
5133  "for_internal_boundary()",
5134  OOMPH_EXCEPTION_LOCATION);
5135  }
5136 #endif
5137 
5138  } // if (!done_segment[is])
5139  } // for (is < nsegments)
5140  } // for (old_is < old_nsegments)
5141 
5142  // For those segments not identified set dummy values, the boundary
5143  // coordinates should be corrected at the synchronisation stage
5144 
5145  // loop over the new segments and check if there not identified
5146  // segments
5147  for (unsigned is = 0; is < nsegments; is++)
5148  {
5149  // Was the segment identified
5150  if (!done_segment[is])
5151  {
5152  // Get the first node of the current segment
5153  FiniteElement* first_seg_ele_pt = segment_sorted_ele_pt[is].front();
5154  // Number of nodes
5155  const unsigned nnod = first_seg_ele_pt->nnode();
5156 
5157  Node* first_seg_node_pt = first_seg_ele_pt->node_pt(0);
5158  if (is_inverted[first_seg_ele_pt])
5159  {
5160  first_seg_node_pt = first_seg_ele_pt->node_pt(nnod - 1);
5161  }
5162 
5163  // Get the arclength for the first node
5164  const double segment_first_node_zeta =
5165  sorted_segment_node_arclength[is][0];
5166 
5167  // Get the node coordinates for the first node
5168  Vector<double> first_node_coord(2);
5169  for (unsigned i = 0; i < 2; i++)
5170  {
5171  first_node_coord[i] = first_seg_node_pt->x(i);
5172  }
5173 
5174  // Get the last node of the current segment
5175  FiniteElement* last_seg_ele_pt = segment_sorted_ele_pt[is].back();
5176  Node* last_seg_node_pt = last_seg_ele_pt->node_pt(nnod - 1);
5177  if (is_inverted[last_seg_ele_pt])
5178  {
5179  last_seg_node_pt = last_seg_ele_pt->node_pt(0);
5180  }
5181 
5182  // Get the arclength for the last node
5183  const double segment_final_node_zeta = segment_arclength[is];
5184 
5185  // Get the node coordinates for the last node
5186  Vector<double> last_node_coord(2);
5187  for (unsigned i = 0; i < 2; i++)
5188  {
5189  last_node_coord[i] = last_seg_node_pt->x(i);
5190  }
5191 
5192  // Re-assign the initial node coordinates
5193  Boundary_segment_initial_coordinate[b][is] = first_node_coord;
5194 
5195  // Check if the boundary has an associated GeomObject
5196  if (this->boundary_geom_object_pt(b) != 0)
5197  {
5198  // Assign the zeta values if the current segment has the
5199  // nodes of the old one
5200 
5201  // If we are in the same order then pass the values as
5202  // they are
5203  Boundary_segment_initial_zeta[b][is] = initial_zeta_segment[is];
5204 
5205  } // if (this->boundary_geom_object_pt(b)!=0)
5206  else
5207  {
5208  // Re-assign the initial arclength for the current segment
5209  Boundary_segment_initial_arclength[b][is] = segment_first_node_zeta;
5210 
5211  } // else if (this->boundary_geom_object_pt(b)!=0)
5212 
5213  // Re-assign the initial node coordinates
5214  Boundary_segment_final_coordinate[b][is] = last_node_coord;
5215 
5216  // Check if the boundary has an associated GeomObject
5217  if (this->boundary_geom_object_pt(b) != 0)
5218  {
5219  // Assign the zeta values if the current segment has the
5220  // nodes of the old one
5221 
5222  // If we are in the same order then pass the values as
5223  // they are
5224  Boundary_segment_final_zeta[b][is] = final_zeta_segment[is];
5225 
5226  } // if (this->boundary_geom_object_pt(b)!=0)
5227  else
5228  {
5229  // Re-assign the final arclength for the current segment
5230  Boundary_segment_final_arclength[b][is] = segment_final_node_zeta;
5231 
5232  } // else if (this->boundary_geom_object_pt(b)!=0)
5233 
5234  Boundary_segment_inverted[b][is] = 0;
5235 
5236  // Mark the segment as done
5237  done_segment[is] = true;
5238 
5239  // Increase the number of re-assigned segments
5240  re_assigned_segments++;
5241 
5242  } // if (!done_segment[is])
5243 
5244  } // for (is < nsegments)
5245 
5246 #ifdef PARANOID
5247  // Compare the number of new segments identified with the old segments
5248  if (re_assigned_segments != nsegments)
5249  {
5250  std::stringstream error_message;
5251  error_message << "Working with boundary (" << b
5252  << ").\nThe number of re-assigned "
5253  << "segments (" << re_assigned_segments
5254  << ") is different from the number\nof segments ("
5255  << nsegments << ")\n\n";
5256  throw OomphLibError(
5257  error_message.str(),
5258  "TriangleMesh::re_assign_initial_zeta_values_for_internal_boundary()",
5259  OOMPH_EXCEPTION_LOCATION);
5260  } // if (re_assigned_segments != nsegments)
5261 #endif
5262 
5263  // Clean all the created face elements
5264  for (unsigned i = 0; i < nele; i++)
5265  {
5266  delete face_el_pt[i];
5267  face_el_pt[i] = 0;
5268  }
5269  }
5270 
5271  /// =====================================================================
5272  /// Select face elements from a given boundary. In case the we are
5273  /// dealing with an internal boundary we use a set of criterias to
5274  /// decide which of the two face elements should be used on represent
5275  /// the internal boundary. We return the face elements, halo or
5276  /// haloed on this processor that form the boundary. The caller method
5277  /// should be in charge of selecting nonhalo elements and deleting the face
5278  /// elements created by this method
5279  /// =====================================================================
5280  template<class ELEMENT>
5282  Vector<FiniteElement*>& face_ele_pt,
5283  const unsigned& b,
5284  bool& is_internal_boundary,
5285  std::map<FiniteElement*, FiniteElement*>& face_to_bulk_element_pt)
5286  {
5287  // Get the communicator of the mesh
5288  OomphCommunicator* comm_pt = this->communicator_pt();
5289 
5290  const unsigned my_rank = comm_pt->my_rank();
5291 
5292  // ------------------------------------------------------------------
5293  // 1) Get the face elements associated with the current boundary
5294  // ------------------------------------------------------------------
5295 
5296  // Temporary storage for face elements (do not take care of
5297  // repeated face elements)
5298  Vector<FiniteElement*> tmp_face_ele_pt;
5299 
5300  const unsigned nregions = this->nregion();
5301 
5302  // If there is more than one region then only use boundary
5303  // coordinates from the bulk side (region 0)
5304  if (nregions > 1)
5305  {
5306  for (unsigned ir = 0; ir < nregions; ir++)
5307  {
5308  const unsigned region_id =
5309  static_cast<unsigned>(this->Region_attribute[ir]);
5310 
5311  // Loop over all elements on boundaries in region -ir-
5312  const unsigned nele_in_region =
5313  this->nboundary_element_in_region(b, region_id);
5314 
5315  // Only bother to do anything else, if there are elements
5316  // associated with the boundary and the current region
5317  if (nele_in_region > 0)
5318  {
5319  // Loop over the bulk elements adjacent to boundary b
5320  for (unsigned e = 0; e < nele_in_region; e++)
5321  {
5322  // Get pointer to the bulk element that is adjacent
5323  // to boundary b
5324  FiniteElement* bulk_ele_pt =
5325  this->boundary_element_in_region_pt(b, region_id, e);
5326 
5327  // Get the index of the face of element e along
5328  // boundary b
5329  int face_index =
5330  this->face_index_at_boundary_in_region(b, region_id, e);
5331 
5332  // Create the face element
5333  FiniteElement* tmp_face_el_pt =
5334  new DummyFaceElement<ELEMENT>(bulk_ele_pt, face_index);
5335 
5336  // Associated the face element with the bulk
5337  face_to_bulk_element_pt[tmp_face_el_pt] = bulk_ele_pt;
5338 
5339  // ... and add it to the tmp storage for all the
5340  // face elements, do not take care for repeated
5341  // ones (at the moment)
5342  tmp_face_ele_pt.push_back(tmp_face_el_pt);
5343 
5344  } // for (e < nele_in_region)
5345 
5346  } // if (nele_in_region > 0)
5347 
5348  } // for (ir < n_regions)
5349 
5350  } // if (n_regions > 1)
5351 
5352  // Otherwise it's just the normal boundary functions
5353  else
5354  {
5355  // Loop over all elements on boundaries
5356  const unsigned nbound_ele = this->nboundary_element(b);
5357 
5358  // Only bother to do anything else, if there are elements
5359  if (nbound_ele > 0)
5360  {
5361  // Loop over the bulk elements adjacent to boundary b
5362  for (unsigned e = 0; e < nbound_ele; e++)
5363  {
5364  // Get pointer to the bulk element that is adjacent to
5365  // boundary b
5366  FiniteElement* bulk_ele_pt = this->boundary_element_pt(b, e);
5367 
5368  // Get the index of the face of element e along
5369  // boundary b
5370  int face_index = this->face_index_at_boundary(b, e);
5371 
5372  // Create the face element
5373  FiniteElement* tmp_face_el_pt =
5374  new DummyFaceElement<ELEMENT>(bulk_ele_pt, face_index);
5375 
5376  // Associated the face element with the bulk
5377  face_to_bulk_element_pt[tmp_face_el_pt] = bulk_ele_pt;
5378 
5379  // ... and add it to the tmp storage for all the face
5380  // elements, do not care for repeated ones (at the
5381  // moment)
5382  tmp_face_ele_pt.push_back(tmp_face_el_pt);
5383 
5384  } // (e < nbound_ele)
5385 
5386  } // (nbound_ele > 0)
5387 
5388  } // else (n_regions > 1)
5389 
5390  // map to know which face element has been already done
5391  std::map<FiniteElement*, bool> done_face;
5392 
5393  // Set the flag to indicate if we are working with an internal
5394  // boundary
5395  is_internal_boundary = false;
5396 
5397  // Free the memory of the elements in this container (only used
5398  // when working with internal boundaries)
5399  Vector<FiniteElement*> free_memory_face_ele_pt;
5400 
5401  // Get the number of face elements in the boundary (including
5402  // repeated)
5403  const unsigned n_tmp_face_ele = tmp_face_ele_pt.size();
5404  for (unsigned ie = 0; ie < n_tmp_face_ele; ie++)
5405  {
5406  // Get the possible main element
5407  FiniteElement* main_face_ele_pt = tmp_face_ele_pt[ie];
5408  if (!done_face[main_face_ele_pt])
5409  {
5410  // Mark the face element as done
5411  done_face[main_face_ele_pt] = true;
5412  // Get the number of nodes for the face element
5413  const unsigned nnodes = main_face_ele_pt->nnode();
5414  // Get the first and last node of the main face element
5415  Node* main_first_node_pt = main_face_ele_pt->node_pt(0);
5416  Node* main_last_node_pt = main_face_ele_pt->node_pt(nnodes - 1);
5417  // Look for the other side face element (we can start from
5418  // the next one, all previous face elements have been
5419  // already identified with its other side face)
5420  for (unsigned iie = ie + 1; iie < n_tmp_face_ele; iie++)
5421  {
5422  // Get the possible dependant element
5423  FiniteElement* dependant_face_ele_pt = tmp_face_ele_pt[iie];
5424  if (!done_face[dependant_face_ele_pt])
5425  {
5426  // Get the first and last node of the dependant
5427  // face element
5428  Node* dependant_first_node_pt = dependant_face_ele_pt->node_pt(0);
5429  Node* dependant_last_node_pt =
5430  dependant_face_ele_pt->node_pt(nnodes - 1);
5431  // Check if the nodes at the ends of both face
5432  // elements match (also check the reversed case)
5433  if (((dependant_first_node_pt == main_first_node_pt) &&
5434  (dependant_last_node_pt == main_last_node_pt)) ||
5435  ((dependant_first_node_pt == main_last_node_pt) &&
5436  (dependant_last_node_pt == main_first_node_pt)))
5437  {
5438  // Set the flag to indicate we are working with an
5439  // internal boundary
5440  is_internal_boundary = true;
5441  // Mark the face element as done
5442  done_face[dependant_face_ele_pt] = true;
5443 
5444  // Now choose which face element will be used
5445  // as the main element. We get the processor in
5446  // charge of the element and choose the one
5447  // with the highest processor in charge or the
5448  // bottom-left bulk element in case the both
5449  // faces are on the same processor
5450 
5451  // Get the bulk element for each face element
5452  // (the main and the dependant face element)
5453  FiniteElement* main_bulk_ele_pt =
5454  face_to_bulk_element_pt[main_face_ele_pt];
5455  FiniteElement* dependant_bulk_ele_pt =
5456  face_to_bulk_element_pt[dependant_face_ele_pt];
5457 
5458  // Get the processor in charge for each bulk
5459  // element
5460  int processor_in_charge_main_bulk_ele =
5461  main_bulk_ele_pt->non_halo_proc_ID();
5462  int processor_in_charge_dependant_bulk_ele =
5463  dependant_bulk_ele_pt->non_halo_proc_ID();
5464 
5465  // If the processor in charge is negative the
5466  // element is not halo, therefore the processor
5467  // in charge is the current one
5468  if (processor_in_charge_main_bulk_ele < 0)
5469  {
5470  processor_in_charge_main_bulk_ele = static_cast<int>(my_rank);
5471  }
5472  if (processor_in_charge_dependant_bulk_ele < 0)
5473  {
5474  processor_in_charge_dependant_bulk_ele =
5475  static_cast<int>(my_rank);
5476  }
5477 
5478  // Flag to know if add the main or dependant
5479  // face element
5480  bool add_main_face_element = true;
5481  if (processor_in_charge_dependant_bulk_ele >
5482  processor_in_charge_main_bulk_ele)
5483  {
5484  // Include the dependant element
5485  add_main_face_element = false;
5486  }
5487  else if (processor_in_charge_main_bulk_ele ==
5488  processor_in_charge_dependant_bulk_ele)
5489  {
5490  // When the processor in charge for both
5491  // elements is the same then use the
5492  // bottom-left criteria on the bulk
5493  // elements to choose the main face element
5494  Vector<double> main_ele_coordinates(2);
5495  Vector<double> dependant_ele_coordinates(2);
5496  // Get the number of nodes on the bulk
5497  // elements
5498  const unsigned n_bulk_nodes = main_bulk_ele_pt->nnode();
5499  for (unsigned inode = 0; inode < n_bulk_nodes; inode++)
5500  {
5501  for (unsigned idim = 0; idim < 2; idim++)
5502  {
5503  main_ele_coordinates[idim] +=
5504  main_bulk_ele_pt->node_pt(inode)->x(idim);
5505  dependant_ele_coordinates[idim] +=
5506  dependant_bulk_ele_pt->node_pt(inode)->x(idim);
5507  } // (idim < 2)
5508 
5509  } // (inode < n_bulk_nodes)
5510 
5511  // Get the average of the nodes coordinates
5512  for (unsigned idim = 0; idim < 2; idim++)
5513  {
5514  main_ele_coordinates[idim] /= (double)n_bulk_nodes;
5515  dependant_ele_coordinates[idim] /= (double)n_bulk_nodes;
5516  }
5517 
5518  // Once we know the average coordinates for
5519  // each element then we choose the one with
5520  // the bottom-left averaged coordinates
5521  if (dependant_ele_coordinates[1] < main_ele_coordinates[1])
5522  {
5523  add_main_face_element = false;
5524  }
5525  else if (dependant_ele_coordinates[1] ==
5526  main_ele_coordinates[1])
5527  {
5528  // The left-most element
5529  if (dependant_ele_coordinates[0] < main_ele_coordinates[0])
5530  {
5531  add_main_face_element = false;
5532  }
5533  }
5534  } // else -- The processor in charge is the
5535  // same for both elements
5536 
5537  if (add_main_face_element)
5538  {
5539  // Add the main face element to the storage
5540  // so we get the halo and haloed nodes from
5541  // it
5542  face_ele_pt.push_back(main_face_ele_pt);
5543  // Mark the dependat face element to free
5544  // its memory
5545  free_memory_face_ele_pt.push_back(dependant_face_ele_pt);
5546  }
5547  else
5548  {
5549  // Add the dependant face element to the
5550  // storage so we get the halo and haloed
5551  // nodes from it
5552  face_ele_pt.push_back(dependant_face_ele_pt);
5553  // Mark the main face element to free its
5554  // memory
5555  free_memory_face_ele_pt.push_back(main_face_ele_pt);
5556  }
5557 
5558  // Break the for to look for the next face
5559  // element
5560  break;
5561 
5562  } // if -- matching of nodes from main ele and
5563  // dependant ele
5564 
5565  } // if (!done_face[dependant_face_ele_pt])
5566 
5567  } // for (iie < n_tmp_face_ele)
5568 
5569  } // if (!done_face[main_face_ele_pt])
5570 
5571  } // for (ie < n_tmp_face_ele)
5572 
5573  // Are there any face element to free its memory
5574  const unsigned n_free_face_ele = free_memory_face_ele_pt.size();
5575  if (n_free_face_ele == 0)
5576  {
5577  // If there is not face elements to free memory that means that
5578  // we are not working with an internal boundary, therefore copy
5579  // all the element from the tmp face elements into the face
5580  // elements container
5581 
5582  // Resize the container
5583  face_ele_pt.resize(n_tmp_face_ele);
5584  // loop over the elements and copy them
5585  for (unsigned i = 0; i < n_tmp_face_ele; i++)
5586  {
5587  face_ele_pt[i] = tmp_face_ele_pt[i];
5588  } // for (i < n_tmp_face_ele)
5589 
5590  } // if (n_free_face_ele == 0)
5591  else
5592  {
5593  // ... otherwise free the memory of the indicated elements
5594  // loop over the elements to free its memory
5595  for (unsigned i = 0; i < n_free_face_ele; i++)
5596  {
5597  delete free_memory_face_ele_pt[i];
5598  free_memory_face_ele_pt[i] = 0;
5599  } // for (i < n_free_face_ele)
5600  }
5601  }
5602 
5603  /// ========================================================================
5604  /// In charge of sinchronize the boundary coordinates for internal
5605  /// boundaries that were split as part of the distribution
5606  /// process. Called after setup_boundary_coordinates() for the
5607  /// original mesh only
5608  /// ========================================================================
5609  template<class ELEMENT>
5611  const unsigned& b)
5612  {
5613  // ------------------------------------------------------------------
5614  // First: Get the face elements associated with the current boundary
5615  // ------------------------------------------------------------------
5616 
5617  // Get the communicator of the mesh
5618  OomphCommunicator* comm_pt = this->communicator_pt();
5619 
5620  const unsigned nproc = comm_pt->nproc();
5621  const unsigned my_rank = comm_pt->my_rank();
5622 
5623  // Temporary storage for face elements (do not take care of repeated
5624  // face elements)
5625  Vector<FiniteElement*> tmp_face_ele_pt;
5626 
5627  const unsigned nregions = this->nregion();
5628 
5629  // map to associate the face element to the bulk element, necessary
5630  // to get the processor in charge for the halo elements
5631  std::map<FiniteElement*, FiniteElement*> face_to_bulk_element_pt;
5632 
5633  // If there is more than one region then only use boundary
5634  // coordinates from the bulk side (region 0)
5635  if (nregions > 1)
5636  {
5637  for (unsigned ir = 0; ir < nregions; ir++)
5638  {
5639  const unsigned region_id =
5640  static_cast<unsigned>(this->Region_attribute[ir]);
5641 
5642  // Loop over all elements on boundaries in region -ir-
5643  const unsigned nele_in_region =
5644  this->nboundary_element_in_region(b, region_id);
5645 
5646  // Only bother to do anything else, if there are elements
5647  // associated with the boundary and the current region
5648  if (nele_in_region > 0)
5649  {
5650  // Loop over the bulk elements adjacent to boundary b
5651  for (unsigned e = 0; e < nele_in_region; e++)
5652  {
5653  // Get pointer to the bulk element that is adjacent to boundary b
5654  FiniteElement* bulk_ele_pt =
5655  this->boundary_element_in_region_pt(b, region_id, e);
5656 
5657  // Get the index of the face of element e along boundary b
5658  int face_index =
5659  this->face_index_at_boundary_in_region(b, region_id, e);
5660 
5661  // Create the face element
5662  FiniteElement* tmp_face_el_pt =
5663  new DummyFaceElement<ELEMENT>(bulk_ele_pt, face_index);
5664 
5665  // ... and add it to the tmp storage for all the face
5666  // elements, do not take care for repeated ones (at the
5667  // moment)
5668  tmp_face_ele_pt.push_back(tmp_face_el_pt);
5669  // Create the map to know if the element is halo
5670  face_to_bulk_element_pt[tmp_face_el_pt] = bulk_ele_pt;
5671 
5672  } // for (e < nele_in_region)
5673 
5674  } // if (nele_in_region > 0)
5675 
5676  } // for (ir < n_regions)
5677 
5678  } // if (n_regions > 1)
5679 
5680  // Otherwise it's just the normal boundary functions
5681  else
5682  {
5683  // Loop over all elements on boundaries
5684  const unsigned nbound_ele = this->nboundary_element(b);
5685 
5686  // Only bother to do anything else, if there are elements
5687  if (nbound_ele > 0)
5688  {
5689  // Loop over the bulk elements adjacent to boundary b
5690  for (unsigned e = 0; e < nbound_ele; e++)
5691  {
5692  // Get pointer to the bulk element that is adjacent to boundary b
5693  FiniteElement* bulk_ele_pt = this->boundary_element_pt(b, e);
5694 
5695  // Get the index of the face of element e along boundary b
5696  int face_index = this->face_index_at_boundary(b, e);
5697 
5698  // Create the face element
5699  FiniteElement* tmp_face_el_pt =
5700  new DummyFaceElement<ELEMENT>(bulk_ele_pt, face_index);
5701 
5702  // ... and add it to the tmp storage for all the face
5703  // elements, do not care for repeated ones (at the moment)
5704  tmp_face_ele_pt.push_back(tmp_face_el_pt);
5705  // Create the map to know if the element is halo
5706  face_to_bulk_element_pt[tmp_face_el_pt] = bulk_ele_pt;
5707 
5708  } // (e < nbound_ele)
5709 
5710  } // (nbound_ele > 0)
5711 
5712  } // else (n_regions > 1)
5713 
5714  // Temporary storage for one side face elements. In case we are
5715  // working with an internal boundary here we store only one of the
5716  // face elements that are at each side of the boundary
5717  Vector<FiniteElement*> face_ele_pt;
5718 
5719  // map to know which face element has been already done
5720  std::map<FiniteElement*, bool> done_face;
5721 
5722  // Flag to indicate if we are working with an internal boundary
5723  bool is_internal_boundary = false;
5724 
5725 #ifdef PARANOID
5726  // Flag to indicate if we are working with an internal boundary (paranoid)
5727  bool is_internal_boundary_paranoid = false;
5728 
5729  // Count the number of other side face elements found in case we are
5730  // working with an internal boundary
5731  unsigned nfound_face_elements = 0;
5732 #endif
5733 
5734  // Get the number of face elements in the boundary
5735  const unsigned nbound_ele = tmp_face_ele_pt.size();
5736  for (unsigned ie = 0; ie < nbound_ele; ie++)
5737  {
5738  // Get the possible main element
5739  FiniteElement* main_face_ele_pt = tmp_face_ele_pt[ie];
5740  if (!done_face[main_face_ele_pt])
5741  {
5742  // Mark the face element as done
5743  done_face[main_face_ele_pt] = true;
5744  // Get the number of nodes for the face element
5745  const unsigned nnodes = main_face_ele_pt->nnode();
5746  // Get the first and last node of the main face element
5747  Node* main_first_node_pt = main_face_ele_pt->node_pt(0);
5748  Node* main_last_node_pt = main_face_ele_pt->node_pt(nnodes - 1);
5749  // Look for the other side face element
5750  for (unsigned iie = ie + 1; iie < nbound_ele; iie++)
5751  {
5752  // Get the possible dependant element
5753  FiniteElement* dependant_face_ele_pt = tmp_face_ele_pt[iie];
5754  if (!done_face[dependant_face_ele_pt])
5755  {
5756  // Get the first and last node of the dependant face element
5757  Node* dependant_first_node_pt = dependant_face_ele_pt->node_pt(0);
5758  Node* dependant_last_node_pt =
5759  dependant_face_ele_pt->node_pt(nnodes - 1);
5760  // Check if the nodes at the ends of both face elements
5761  // match (also check the reversed case)
5762  if (((dependant_first_node_pt == main_first_node_pt) &&
5763  (dependant_last_node_pt == main_last_node_pt)) ||
5764  ((dependant_first_node_pt == main_last_node_pt) &&
5765  (dependant_last_node_pt == main_first_node_pt)))
5766  {
5767 #ifdef PARANOID
5768  // Increase the number of found face elements
5769  nfound_face_elements += 2;
5770 #endif
5771  // Set the flag to indicate we are working with an
5772  // internal boundary
5773  is_internal_boundary = true;
5774  // Mark the face element as done
5775  done_face[dependant_face_ele_pt] = true;
5776 
5777  // Now choose which face element will be used as the main
5778  // element. Use the same criteria as the compute segments
5779  // connectivity method (highest processor in charge or
5780  // bottom-left bulk element)
5781 
5782  // Get the bulk element for each face element (the main
5783  // and the dependant face element)
5784  FiniteElement* main_bulk_ele_pt =
5785  face_to_bulk_element_pt[main_face_ele_pt];
5786  FiniteElement* dependant_bulk_ele_pt =
5787  face_to_bulk_element_pt[dependant_face_ele_pt];
5788 
5789  // Get the processor in charge for each bulk element
5790  int processor_in_charge_main_bulk_ele =
5791  main_bulk_ele_pt->non_halo_proc_ID();
5792  int processor_in_charge_dependant_bulk_ele =
5793  dependant_bulk_ele_pt->non_halo_proc_ID();
5794 
5795  // If the processor in charge is negative the element is
5796  // not halo, therefore the processor in charge is the
5797  // current one
5798  if (processor_in_charge_main_bulk_ele < 0)
5799  {
5800  processor_in_charge_main_bulk_ele = static_cast<int>(my_rank);
5801  }
5802  if (processor_in_charge_dependant_bulk_ele < 0)
5803  {
5804  processor_in_charge_dependant_bulk_ele =
5805  static_cast<int>(my_rank);
5806  }
5807 
5808  // Flag to know if add the main or dependant face element
5809  bool add_main_face_element = true;
5810  if (processor_in_charge_dependant_bulk_ele >
5811  processor_in_charge_main_bulk_ele)
5812  {
5813  // Include the dependant element
5814  add_main_face_element = false;
5815  }
5816  else if (processor_in_charge_main_bulk_ele ==
5817  processor_in_charge_dependant_bulk_ele)
5818  {
5819  // When the processor in charge for both elements is the same
5820  // then use the bottom-left criteria on the bulk elements to
5821  // choose the main face element
5822  Vector<double> main_ele_coordinates(2);
5823  Vector<double> dependant_ele_coordinates(2);
5824  // Get the number of nodes on the bulk elements
5825  const unsigned n_bulk_nodes = main_bulk_ele_pt->nnode();
5826  for (unsigned inode = 0; inode < n_bulk_nodes; inode++)
5827  {
5828  for (unsigned idim = 0; idim < 2; idim++)
5829  {
5830  main_ele_coordinates[idim] +=
5831  main_bulk_ele_pt->node_pt(inode)->x(idim);
5832  dependant_ele_coordinates[idim] +=
5833  dependant_bulk_ele_pt->node_pt(inode)->x(idim);
5834  } // (idim < 2)
5835  } // (inode < n_bulk_nodes)
5836 
5837  // Get the average of the nodes coordinates
5838  for (unsigned idim = 0; idim < 2; idim++)
5839  {
5840  main_ele_coordinates[idim] /= (double)n_bulk_nodes;
5841  dependant_ele_coordinates[idim] /= (double)n_bulk_nodes;
5842  }
5843 
5844  // Once we know the average coordinates for each element
5845  // then we choose the one with the bottom-left averaged
5846  // coordinates
5847  if (dependant_ele_coordinates[1] < main_ele_coordinates[1])
5848  {
5849  add_main_face_element = false;
5850  }
5851  else if (dependant_ele_coordinates[1] ==
5852  main_ele_coordinates[1])
5853  {
5854  // The left-most element
5855  if (dependant_ele_coordinates[0] < main_ele_coordinates[0])
5856  {
5857  add_main_face_element = false;
5858  }
5859  }
5860  } // else -- The processor in charge is the same for both
5861  // elements
5862 
5863  if (add_main_face_element)
5864  {
5865  // Add the main face element to the storage so we get
5866  // the halo and haloed nodes from these face element
5867  face_ele_pt.push_back(main_face_ele_pt);
5868  }
5869  else
5870  {
5871  // Add the main face element to the storage so we get
5872  // the halo and haloed nodes from these face element
5873  face_ele_pt.push_back(dependant_face_ele_pt);
5874  }
5875 
5876  // Break the for to look for the next face element
5877  break;
5878 
5879  } // if -- matching of nodes from main ele and dependant ele
5880  } // if (!done_face[dependant_face_ele_pt])
5881  } // for (iie < nbound_ele)
5882  } // if (!done_face[main_face_ele_pt])
5883  } // for (ie < nbound_ele)
5884 
5885  // Get the number of face elements
5886  const unsigned nface_ele = face_ele_pt.size();
5887 
5888 #ifdef PARANOID
5889  // Check if we are working with an internal open curve. First check
5890  // if there are elements, in a distributed approach they may be no
5891  // elements associated to the boundary
5892  if (nbound_ele > 0 && nfound_face_elements == nbound_ele)
5893  {
5894  is_internal_boundary_paranoid = true;
5895  }
5896 
5897  if (nbound_ele > 0 && is_internal_boundary_paranoid &&
5898  nbound_ele != nface_ele * 2)
5899  {
5900  std::ostringstream error_message;
5901  error_message
5902  << "The info. to perform the synchronisation of the boundary "
5903  << "coordinates was not completely established\n"
5904  << "In this case it was the number of non repeated boundary elements\n"
5905  << "Number of boundary elements: (" << nbound_ele << ")\n"
5906  << "Number of nonrepeated boundary elements: (" << nface_ele << ")\n";
5907  throw OomphLibError(error_message.str(),
5908  "TriangleMesh::synchronize_boundary_coordinates()",
5909  OOMPH_EXCEPTION_LOCATION);
5910  }
5911 #endif
5912 
5913  // ----------------------------------------------------------------
5914  // Second: Identify the halo face elements
5915  // ----------------------------------------------------------------
5916 
5917  // A flag vector to mark those face elements that are considered as
5918  // halo in the current processor
5919  std::vector<bool> is_halo_face_element(nface_ele, false);
5920 
5921  // Count the total number of non halo face elements
5922  unsigned nnon_halo_face_elements = 0;
5923 
5924  for (unsigned ie = 0; ie < nface_ele; ie++)
5925  {
5926  FiniteElement* face_el_pt = face_ele_pt[ie];
5927  // Get the bulk element
5928  FiniteElement* tmp_bulk_ele_pt = face_to_bulk_element_pt[face_el_pt];
5929  // Check if the bulk element is halo
5930  if (!tmp_bulk_ele_pt->is_halo())
5931  {
5932  is_halo_face_element[ie] = false;
5933  nnon_halo_face_elements++;
5934  }
5935  else
5936  {
5937  // Mark the face element as halo
5938  is_halo_face_element[ie] = true;
5939  }
5940  } // for (ie < nface_ele)
5941 
5942  // -----------------------------------------------------------------
5943  // Third: Go through the face elements and get the nodes from the
5944  // elements. The boundary coordinate from each node is sent to its
5945  // processor in charge, then that processor will be responsible to
5946  // send the bound coordinate to all the processors that have a halo
5947  // representation of the node
5948  // -----------------------------------------------------------------
5949 
5950  // A map to know which nodes are already done
5951  std::map<Node*, bool> done_node;
5952 
5953  // The storage for the halo nodes on face elements in this processor
5954  // with other processors
5955  Vector<Vector<Node*>> face_halo_node_pt(nproc);
5956 
5957  // The storage for the ids of the halo nodes on face elements in
5958  // this processor with other processors
5959  Vector<Vector<unsigned>> face_halo_node_id(nproc);
5960 
5961  // The storage for the haloed nodes on face elements in this
5962  // processor with other processors
5963  Vector<Vector<Node*>> face_haloed_node_pt(nproc);
5964 
5965  // The storage for the ids of the haloed nodes on face elements in
5966  // this processor with other processors
5967  Vector<Vector<unsigned>> face_haloed_node_id(nproc);
5968 
5969  // A map to know which nodes are face nodes and the processor in
5970  // charge is the current one
5971  std::map<Node*, bool> done_haloed_face_node;
5972 
5973  // Go through all the face elements
5974  for (unsigned iface = 0; iface < nface_ele; iface++)
5975  {
5976  // Only work with the non halo face elements
5977  if (!is_halo_face_element[iface])
5978  {
5979  // Get the face element
5980  FiniteElement* ele_face_pt = face_ele_pt[iface];
5981  // The number of nodes of the face elements
5982  const unsigned nnodes = ele_face_pt->nnode();
5983  // Go through all the nodes in the face element
5984  for (unsigned in = 0; in < nnodes; in++)
5985  {
5986  Node* face_node_pt = ele_face_pt->node_pt(in);
5987  // Check if node is done
5988  if (!done_node[face_node_pt])
5989  {
5990  // Mark the node as done
5991  done_node[face_node_pt] = true;
5992  // First check if the node is halo
5993  if (face_node_pt->is_halo())
5994  {
5995  // Get the processor in charge for the current node
5996  int int_nonhalo_ID = face_node_pt->non_halo_proc_ID();
5997 #ifdef PARANOID
5998  if (int_nonhalo_ID < 0)
5999  {
6000  std::ostringstream error_message;
6001  error_message
6002  << "The node was marked to be halo but the processor in "
6003  << "charge was found to be -1\n\n";
6004  throw OomphLibError(
6005  error_message.str(),
6006  "TriangleMesh::synchronize_boundary_coordinates()",
6007  OOMPH_EXCEPTION_LOCATION);
6008  }
6009 #endif
6010  const unsigned ip = static_cast<unsigned>(int_nonhalo_ID);
6011  // Add the node to the structure that holds the halo
6012  // nodes, the current processor will need to send the
6013  // info. to the processor in charge.
6014  face_halo_node_pt[ip].push_back(face_node_pt);
6015  // ... finally look for the halo id with the processor in
6016  // charge
6017 #ifdef PARANOID
6018  bool found_halo_node = false;
6019 #endif
6020  const unsigned nhalo_iproc = this->nhalo_node(ip);
6021  for (unsigned ihn = 0; ihn < nhalo_iproc; ihn++)
6022  {
6023  Node* compare_face_node_pt = this->halo_node_pt(ip, ihn);
6024  if (compare_face_node_pt == face_node_pt)
6025  {
6026  // Once found the id of the node with the processor
6027  // store the id in the proper storage
6028  face_halo_node_id[ip].push_back(ihn);
6029 #ifdef PARANOID
6030  // Set the flag to mark as found the halo node
6031  found_halo_node = true;
6032 #endif
6033  // Break the loop
6034  break;
6035  }
6036  } // for (ih < nhalo_iproc)
6037 #ifdef PARANOID
6038  if (!found_halo_node)
6039  {
6040  std::ostringstream error_message;
6041  error_message
6042  << "The halo id of the current node: (" << face_node_pt->x(0)
6043  << ", " << face_node_pt->x(1) << ") with processor (" << ip
6044  << ") was not found!!!\n\n";
6045  throw OomphLibError(
6046  error_message.str(),
6047  "TriangleMesh::synchronize_boundary_coordinates()",
6048  OOMPH_EXCEPTION_LOCATION);
6049  }
6050 #endif
6051  } // if (face_node_pt->is_halo())
6052  // If the node is not halo then it could be haloed. If that
6053  // is the case then store the processors at which the node
6054  // is haloed and its id. The info. of these nodes will be
6055  // sent to all the processors with a halo counterpart
6056  else
6057  {
6058  for (unsigned ip = 0; ip < nproc; ip++)
6059  {
6060  // Only work with processors different that the current one
6061  if (ip != my_rank)
6062  {
6063  // If the node is found to be haloed with the "ip"
6064  // processor then save the haloed id in the storage.
6065  // The current processor needs to send info. to the
6066  // other processors to establish the boundary
6067  // coordinates
6068 
6069  // Get the number of haloed nodes with processor ip
6070  const unsigned nhaloed_iproc = this->nhaloed_node(ip);
6071  for (unsigned ihdn = 0; ihdn < nhaloed_iproc; ihdn++)
6072  {
6073  Node* compare_face_node_pt = this->haloed_node_pt(ip, ihdn);
6074  if (face_node_pt == compare_face_node_pt)
6075  {
6076  // Store the node on the haloed node vector for
6077  // the corresponding processor
6078  face_haloed_node_pt[ip].push_back(face_node_pt);
6079  // Now store the halo id of the node with the
6080  // current processor
6081  face_haloed_node_id[ip].push_back(ihdn);
6082  // Mark the node as haloed with other processors,
6083  // so we know the processor in charge is the
6084  // current one "my_rank".
6085  done_haloed_face_node[face_node_pt] = true;
6086  // Break looking in the current processor, look in
6087  // the next one
6088  break;
6089  } // if (face_node_pt == compare_face_node_pt)
6090  } // for (ihdn < nhaloed_node_iproc)
6091  } // if (ip != my_rank)
6092  } // for (ip < nproc)
6093  } // else (non halo node)
6094  } // if (!done_node[node_face_pt])
6095  } // for (in < nnodes)
6096  } // if (!is_halo_face_element[iface])
6097  } // for (iface < nface_ele)
6098 
6099  // -----------------------------------------------------------------
6100  // Fourth: Go through the halo nodes, package and send the
6101  // info. necessary to identify the face nodes in the processor in
6102  // charge. Identify the haloed nodes in the processor in charge and
6103  // establish the boundary coordinates, check if those nodes are
6104  // (already) marked as faced nodes, if that is the case then do not
6105  // establish the boundary coordinates but register them to send back
6106  // the info. to all the processors that have a halo representation
6107  // of the face node
6108  // -----------------------------------------------------------------
6109 
6110  // Go through all processors
6111  for (unsigned ip = 0; ip < nproc; ip++)
6112  {
6113  // Only work with processors different than the current one
6114  if (ip != my_rank)
6115  {
6116  const unsigned nhalo_face_nodes = face_halo_node_pt[ip].size();
6117 #ifdef PARANOID
6118  if (nhalo_face_nodes != face_halo_node_id[ip].size())
6119  {
6120  std::ostringstream error_message;
6121  error_message
6122  << "The number of found halo face nodes (" << nhalo_face_nodes
6123  << ") is different from the number of\nfound halo face ids ("
6124  << face_halo_node_id[ip].size() << ")!!!\n\n";
6125  throw OomphLibError(
6126  error_message.str(),
6127  "TriangleMesh::synchronize_boundary_coordinates()",
6128  OOMPH_EXCEPTION_LOCATION);
6129  }
6130 #endif
6131 
6132  // Container to send the info. related with the halo nodes to be
6133  // identified in the processors in charge
6134  Vector<unsigned> flat_unsigned_send_packed_data;
6135  Vector<double> flat_double_send_packed_data;
6136 
6137  // Go through the halo face nodes in the "ip" processor
6138  for (unsigned ihfn = 0; ihfn < nhalo_face_nodes; ihfn++)
6139  {
6140  // Get the "ihfn"-th face node with the "ip" processor
6141  Node* halo_face_node_pt = face_halo_node_pt[ip][ihfn];
6142  // Get the halo id with the "ip" processor
6143  const unsigned halo_id = face_halo_node_id[ip][ihfn];
6144  // Get the boundary coordinate of the node
6145  Vector<double> zeta(1);
6146  halo_face_node_pt->get_coordinates_on_boundary(b, zeta);
6147  // Store the info. in the containers
6148  flat_unsigned_send_packed_data.push_back(halo_id);
6149  flat_double_send_packed_data.push_back(zeta[0]);
6150  }
6151 
6152  // Send the info.
6153  MPI_Status status;
6154  MPI_Request request;
6155 
6156  // Processor to which send the info
6157  int send_proc = static_cast<int>(ip);
6158  // Processor from which receive the info
6159  int receive_proc = static_cast<int>(ip);
6160 
6161  // Storage to receive the info.
6162  Vector<unsigned> flat_unsigned_receive_packed_data;
6163  Vector<double> flat_double_receive_packed_data;
6164 
6165  // --------------
6166  // Unsigned data
6167  unsigned nflat_unsigned_send = flat_unsigned_send_packed_data.size();
6168  MPI_Isend(&nflat_unsigned_send,
6169  1,
6170  MPI_UNSIGNED,
6171  send_proc,
6172  1,
6173  comm_pt->mpi_comm(),
6174  &request);
6175 
6176  unsigned nflat_unsigned_receive = 0;
6177  MPI_Recv(&nflat_unsigned_receive,
6178  1,
6179  MPI_UNSIGNED,
6180  receive_proc,
6181  1,
6182  comm_pt->mpi_comm(),
6183  &status);
6184 
6185  MPI_Wait(&request, MPI_STATUS_IGNORE);
6186 
6187  if (nflat_unsigned_send != 0)
6188  {
6189  MPI_Isend(&flat_unsigned_send_packed_data[0],
6190  nflat_unsigned_send,
6191  MPI_UNSIGNED,
6192  send_proc,
6193  2,
6194  comm_pt->mpi_comm(),
6195  &request);
6196  }
6197 
6198  if (nflat_unsigned_receive != 0)
6199  {
6200  flat_unsigned_receive_packed_data.resize(nflat_unsigned_receive);
6201  MPI_Recv(&flat_unsigned_receive_packed_data[0],
6202  nflat_unsigned_receive,
6203  MPI_UNSIGNED,
6204  receive_proc,
6205  2,
6206  comm_pt->mpi_comm(),
6207  &status);
6208  }
6209 
6210  if (nflat_unsigned_send != 0)
6211  {
6212  MPI_Wait(&request, MPI_STATUS_IGNORE);
6213  }
6214 
6215  // --------------
6216  // Double data
6217  unsigned nflat_double_send = flat_double_send_packed_data.size();
6218  MPI_Isend(&nflat_double_send,
6219  1,
6220  MPI_DOUBLE,
6221  send_proc,
6222  3,
6223  comm_pt->mpi_comm(),
6224  &request);
6225 
6226  unsigned nflat_double_receive = 0;
6227  MPI_Recv(&nflat_double_receive,
6228  1,
6229  MPI_DOUBLE,
6230  receive_proc,
6231  3,
6232  comm_pt->mpi_comm(),
6233  &status);
6234 
6235  MPI_Wait(&request, MPI_STATUS_IGNORE);
6236 
6237  if (nflat_double_send != 0)
6238  {
6239  MPI_Isend(&flat_double_send_packed_data[0],
6240  nflat_double_send,
6241  MPI_DOUBLE,
6242  send_proc,
6243  4,
6244  comm_pt->mpi_comm(),
6245  &request);
6246  }
6247 
6248  if (nflat_double_receive != 0)
6249  {
6250  flat_double_receive_packed_data.resize(nflat_double_receive);
6251  MPI_Recv(&flat_double_receive_packed_data[0],
6252  nflat_double_receive,
6253  MPI_DOUBLE,
6254  receive_proc,
6255  4,
6256  comm_pt->mpi_comm(),
6257  &status);
6258  }
6259 
6260  if (nflat_double_send != 0)
6261  {
6262  MPI_Wait(&request, MPI_STATUS_IGNORE);
6263  }
6264  // --------------
6265 
6266 #ifdef PARANOID
6267  if (nflat_unsigned_receive != nflat_double_receive)
6268  {
6269  std::ostringstream error_message;
6270  error_message << "The number of unsigned received data ("
6271  << nflat_unsigned_receive << ") is different from the "
6272  << "number\nof double received data ("
6273  << nflat_double_receive << ")!!!\n\n";
6274  throw OomphLibError(
6275  error_message.str(),
6276  "TriangleMesh::synchronize_boundary_coordinates()",
6277  OOMPH_EXCEPTION_LOCATION);
6278  }
6279 #endif
6280 
6281  // With the received info. establish the boundary coordinates
6282  // for the face nodes that this processor is in charge (haloed
6283  // nodes)
6284  for (unsigned iflat_packed = 0; iflat_packed < nflat_unsigned_receive;
6285  iflat_packed++)
6286  {
6287  // Get the haloed id for the node
6288  const unsigned haloed_id =
6289  flat_unsigned_receive_packed_data[iflat_packed];
6290  // Get the boundary coordinates
6291  Vector<double> zeta(1);
6292  zeta[0] = flat_double_receive_packed_data[iflat_packed];
6293 
6294  // Get the haloed node
6295  Node* haloed_face_node_pt = this->haloed_node_pt(ip, haloed_id);
6296 
6297  // If the node has already set the boundary coordinates then
6298  // do not establish it. This is the case for the nodes that
6299  // lie on the boundary, for those nodes not identified on the
6300  // boundary since no elements lie on the boundary but the node
6301  // is on the boundary (a corner of an element lies on the
6302  // boundary) set boundary coordinates and register them to
6303  // send their info. to the processors with a halo counterpart
6304 
6305  // If the node is not haloed face in the procesor in charge
6306  // then set the boundary coordinates and register the node to
6307  // send back the boundary coordinates to the processors with a
6308  // halo counterpart
6309  if (!done_haloed_face_node[haloed_face_node_pt])
6310  {
6311  // Establish the boundary coordinates
6312  haloed_face_node_pt->set_coordinates_on_boundary(b, zeta);
6313 
6314  // Look in all processors where the node could be halo
6315  for (unsigned iiproc = 0; iiproc < nproc; iiproc++)
6316  {
6317  // Only work with processors different than the current one
6318  if (iiproc != my_rank)
6319  {
6320  // Get the number of haloed nodes with processor iiproc
6321  const unsigned nhaloed_node_iiproc = this->nhaloed_node(iiproc);
6322  for (unsigned ihdn = 0; ihdn < nhaloed_node_iiproc; ihdn++)
6323  {
6324  Node* compare_haloed_node_pt =
6325  this->haloed_node_pt(iiproc, ihdn);
6326  if (haloed_face_node_pt == compare_haloed_node_pt)
6327  {
6328  // Store the node on the haloed node vector for the
6329  // corresponding processor
6330  face_haloed_node_pt[iiproc].push_back(haloed_face_node_pt);
6331  // Now store the halo id of the node with the current
6332  // processor
6333  face_haloed_node_id[iiproc].push_back(ihdn);
6334  // Break searching in the current processor, search in
6335  // the next one
6336  break;
6337  } // if (haloed_face_node_pt==compare_haloed_face_node_pt)
6338  } // for (ihdn < nhaloed_node_iproc)
6339  } // if (iiproc != my_rank)
6340  } // for (iiproc < nproc)
6341  } // if (!done_haloed_face_node[haloed_face_node_pt])
6342  } // for (iflat_packed < nflat_unsigned_receive)
6343  } // if (ip != my_rank)
6344  } // for (ip < nproc)
6345 
6346  // -----------------------------------------------------------------
6347  // Fifth: The boundary coordinates have been established in the
6348  // processors in charge of the nodes. Now each processor send back
6349  // the boundary coordinates to all the processors where there is a
6350  // halo representation of the node
6351  // -----------------------------------------------------------------
6352 
6353  // Go through all processors
6354  for (unsigned ip = 0; ip < nproc; ip++)
6355  {
6356  // Only work with processors different than the current one
6357  if (ip != my_rank)
6358  {
6359  // Container to send the info. of the haloed nodes to all the
6360  // processors
6361  Vector<unsigned> flat_unsigned_send_packed_data;
6362  Vector<double> flat_double_send_packed_data;
6363 
6364  // Get the total number of haloed face nodes with the "ip"
6365  // processor
6366  const unsigned nhaloed_face_nodes = face_haloed_node_pt[ip].size();
6367  // Go through the haloed face nodes in the "ip" processor
6368  for (unsigned ihdfn = 0; ihdfn < nhaloed_face_nodes; ihdfn++)
6369  {
6370  // Get the "ihdfn"-th face node with the "ip" processor
6371  Node* haloed_face_node_pt = face_haloed_node_pt[ip][ihdfn];
6372  // Get the haloed id with the "ip" processor
6373  const unsigned haloed_id = face_haloed_node_id[ip][ihdfn];
6374  // Get the boundary coordinate of the node
6375  Vector<double> zeta(1);
6376  haloed_face_node_pt->get_coordinates_on_boundary(b, zeta);
6377  // Store the info. in the containers
6378  flat_unsigned_send_packed_data.push_back(haloed_id);
6379  flat_double_send_packed_data.push_back(zeta[0]);
6380  }
6381 
6382  // Send the info.
6383  MPI_Status status;
6384  MPI_Request request;
6385 
6386  // Processor to which send the info
6387  int send_proc = static_cast<int>(ip);
6388  // Processor from which receive the info
6389  int receive_proc = static_cast<int>(ip);
6390 
6391  // Storage to receive the info.
6392  Vector<unsigned> flat_unsigned_receive_packed_data;
6393  Vector<double> flat_double_receive_packed_data;
6394 
6395  // --------------
6396  // Unsigned data
6397  unsigned nflat_unsigned_send = flat_unsigned_send_packed_data.size();
6398  MPI_Isend(&nflat_unsigned_send,
6399  1,
6400  MPI_UNSIGNED,
6401  send_proc,
6402  1,
6403  comm_pt->mpi_comm(),
6404  &request);
6405 
6406  unsigned nflat_unsigned_receive = 0;
6407  MPI_Recv(&nflat_unsigned_receive,
6408  1,
6409  MPI_UNSIGNED,
6410  receive_proc,
6411  1,
6412  comm_pt->mpi_comm(),
6413  &status);
6414 
6415  MPI_Wait(&request, MPI_STATUS_IGNORE);
6416 
6417  if (nflat_unsigned_send != 0)
6418  {
6419  MPI_Isend(&flat_unsigned_send_packed_data[0],
6420  nflat_unsigned_send,
6421  MPI_UNSIGNED,
6422  send_proc,
6423  2,
6424  comm_pt->mpi_comm(),
6425  &request);
6426  }
6427 
6428  if (nflat_unsigned_receive != 0)
6429  {
6430  flat_unsigned_receive_packed_data.resize(nflat_unsigned_receive);
6431  MPI_Recv(&flat_unsigned_receive_packed_data[0],
6432  nflat_unsigned_receive,
6433  MPI_UNSIGNED,
6434  receive_proc,
6435  2,
6436  comm_pt->mpi_comm(),
6437  &status);
6438  }
6439 
6440  if (nflat_unsigned_send != 0)
6441  {
6442  MPI_Wait(&request, MPI_STATUS_IGNORE);
6443  }
6444 
6445  // --------------
6446  // Double data
6447  unsigned nflat_double_send = flat_double_send_packed_data.size();
6448  MPI_Isend(&nflat_double_send,
6449  1,
6450  MPI_DOUBLE,
6451  send_proc,
6452  3,
6453  comm_pt->mpi_comm(),
6454  &request);
6455 
6456  unsigned nflat_double_receive = 0;
6457  MPI_Recv(&nflat_double_receive,
6458  1,
6459  MPI_DOUBLE,
6460  receive_proc,
6461  3,
6462  comm_pt->mpi_comm(),
6463  &status);
6464 
6465  MPI_Wait(&request, MPI_STATUS_IGNORE);
6466 
6467  if (nflat_double_send != 0)
6468  {
6469  MPI_Isend(&flat_double_send_packed_data[0],
6470  nflat_double_send,
6471  MPI_DOUBLE,
6472  send_proc,
6473  4,
6474  comm_pt->mpi_comm(),
6475  &request);
6476  }
6477 
6478  if (nflat_double_receive != 0)
6479  {
6480  flat_double_receive_packed_data.resize(nflat_double_receive);
6481  MPI_Recv(&flat_double_receive_packed_data[0],
6482  nflat_double_receive,
6483  MPI_DOUBLE,
6484  receive_proc,
6485  4,
6486  comm_pt->mpi_comm(),
6487  &status);
6488  }
6489 
6490  if (nflat_double_send != 0)
6491  {
6492  MPI_Wait(&request, MPI_STATUS_IGNORE);
6493  }
6494  // --------------
6495 
6496 #ifdef PARANOID
6497  if (nflat_unsigned_receive != nflat_double_receive)
6498  {
6499  std::ostringstream error_message;
6500  error_message << "The number of unsigned received data ("
6501  << nflat_unsigned_receive << ") is different from the "
6502  << "number\nof double received data ("
6503  << nflat_double_receive << ")!!!\n\n";
6504  throw OomphLibError(
6505  error_message.str(),
6506  "TriangleMesh::synchronize_boundary_coordinates()",
6507  OOMPH_EXCEPTION_LOCATION);
6508  }
6509 #endif
6510 
6511  // With the received info. establish the boundary coordinates
6512  // received for the face nodes that this processor is not in
6513  // charge (halo nodes)
6514  for (unsigned iflat_packed = 0; iflat_packed < nflat_unsigned_receive;
6515  iflat_packed++)
6516  {
6517  // Get the halo id for the node
6518  const unsigned halo_id =
6519  flat_unsigned_receive_packed_data[iflat_packed];
6520  // Get the boundary coordinates
6521  Vector<double> zeta(1);
6522  zeta[0] = flat_double_receive_packed_data[iflat_packed];
6523 
6524  // Get the halo node
6525  Node* halo_face_node_pt = this->halo_node_pt(ip, halo_id);
6526 
6527  // It could be possible that the node has been already
6528  // established boundary coordinates since it is a halo face
6529  // node. However, for those elements not on the boundary, but
6530  // having a corner node on the boundary this procedure will
6531  // establish boundary coordinates for those nodes
6532 
6533  // this->add_boundary_node(b, halo_face_node_pt);
6534 
6535  // Establish the boundary coordinates
6536  halo_face_node_pt->set_coordinates_on_boundary(b, zeta);
6537  } // for (iflat_packed < nflat_unsigned_receive)
6538  } // if (ip != my_rank)
6539  } // for (ip < nproc)
6540 
6541  // Clean all the created face elements
6542  for (unsigned ie = 0; ie < nbound_ele; ie++)
6543  {
6544  delete tmp_face_ele_pt[ie];
6545  tmp_face_ele_pt[ie] = 0;
6546  }
6547 
6548  // Now get a new face mesh representation and fill the data for those
6549  // processors with halo segments
6550  if (is_internal_boundary)
6551  {
6552  re_scale_re_assigned_initial_zeta_values_for_internal_boundary(b);
6553  }
6554  }
6555 
6556  //======================================================================
6557  /// Re-assign the boundary segments initial zeta (arclength)
6558  /// for those internal boundaries that were splited during the
6559  /// distribution process (only apply for internal boundaries that
6560  /// have one face element at each side of the boundary)
6561  //======================================================================
6562  template<class ELEMENT>
6565  const unsigned& b)
6566  {
6567  // ------------------------------------------------------------------
6568  // First: Get the face elements associated with the current boundary
6569  // Only include nonhalo face elements
6570  // ------------------------------------------------------------------
6571  // Temporary storage for face elements
6572  Vector<FiniteElement*> face_el_pt;
6573 
6574  // Temporary storage for the number of elements adjacent to the
6575  // boundary
6576  unsigned nele = 0;
6577 
6578  // Temporary storage for elements adjacent to the boundary that have
6579  // a common edge (related with internal boundaries)
6580  unsigned n_repeated_ele = 0;
6581 
6582  const unsigned n_regions = this->nregion();
6583 
6584  // Temporary storage for already done nodes
6585  Vector<std::pair<Node*, Node*>> done_nodes_pt;
6586 
6587  // If there is more than one region then only use boundary
6588  // coordinates from the bulk side (region 0)
6589  if (n_regions > 1)
6590  {
6591  for (unsigned rr = 0; rr < n_regions; rr++)
6592  {
6593  const unsigned region_id =
6594  static_cast<unsigned>(this->Region_attribute[rr]);
6595 
6596  // Loop over all elements on boundaries in region i_r
6597  const unsigned nel_in_region =
6598  this->nboundary_element_in_region(b, region_id);
6599 
6600  unsigned nel_repetead_in_region = 0;
6601 
6602  // Only bother to do anything else, if there are elements
6603  // associated with the boundary and the current region
6604  if (nel_in_region > 0)
6605  {
6606  bool repeated = false;
6607 
6608  // Loop over the bulk elements adjacent to boundary b
6609  for (unsigned e = 0; e < nel_in_region; e++)
6610  {
6611  // Get pointer to the bulk element that is adjacent to
6612  // boundary b
6613  FiniteElement* bulk_elem_pt =
6614  this->boundary_element_in_region_pt(b, region_id, e);
6615 
6616  // Remember only to work with nonhalo elements
6617  if (bulk_elem_pt->is_halo())
6618  {
6619  n_repeated_ele++;
6620  continue;
6621  }
6622 
6623  // Find the index of the face of element e along boundary b
6624  int face_index =
6625  this->face_index_at_boundary_in_region(b, region_id, e);
6626 
6627  // Before adding the new element we need to be sure that the
6628  // edge that this element represent has not been already
6629  // added
6630  FiniteElement* tmp_ele_pt =
6631  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
6632 
6633  const unsigned n_nodes = tmp_ele_pt->nnode();
6634 
6635  std::pair<Node*, Node*> tmp_pair = std::make_pair(
6636  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
6637 
6638  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
6639  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
6640 
6641  // Search for repeated nodes
6642  const unsigned n_done_nodes = done_nodes_pt.size();
6643  for (unsigned l = 0; l < n_done_nodes; l++)
6644  {
6645  if (tmp_pair == done_nodes_pt[l] ||
6646  tmp_pair_inverse == done_nodes_pt[l])
6647  {
6648  nel_repetead_in_region++;
6649  repeated = true;
6650  break;
6651  }
6652  }
6653 
6654  // Create new face element
6655  if (!repeated)
6656  {
6657  // Add the pair of nodes (edge) to the node dones
6658  done_nodes_pt.push_back(tmp_pair);
6659  // Add the element to the face elements
6660  face_el_pt.push_back(tmp_ele_pt);
6661  }
6662  else
6663  {
6664  // Clean up
6665  delete tmp_ele_pt;
6666  tmp_ele_pt = 0;
6667  }
6668 
6669  // Re-start
6670  repeated = false;
6671 
6672  } // for nel
6673 
6674  nele += nel_in_region;
6675 
6676  n_repeated_ele += nel_repetead_in_region;
6677 
6678  } // if (nel_in_region > 0)
6679  } // for (rr < n_regions)
6680  } // if (n_regions > 1)
6681  // Otherwise it's just the normal boundary functions
6682  else
6683  {
6684  // Loop over all elements on boundaries
6685  nele = this->nboundary_element(b);
6686 
6687  // Only bother to do anything else, if there are elements
6688  if (nele > 0)
6689  {
6690  // Check for repeated ones
6691  bool repeated = false;
6692 
6693  // Loop over the bulk elements adjacent to boundary b
6694  for (unsigned e = 0; e < nele; e++)
6695  {
6696  // Get pointer to the bulk element that is adjacent to
6697  // boundary b
6698  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
6699 
6700  // Remember only to work with nonhalo elements
6701  if (bulk_elem_pt->is_halo())
6702  {
6703  n_repeated_ele++;
6704  // Skip the halo element
6705  continue;
6706  }
6707 
6708  // Find the index of the face of element e along boundary b
6709  int face_index = this->face_index_at_boundary(b, e);
6710 
6711  // Before adding the new element we need to be sure that the
6712  // edge that this element represents has not been already
6713  // added (only applies for internal boundaries)
6714  FiniteElement* tmp_ele_pt =
6715  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
6716 
6717  const unsigned n_nodes = tmp_ele_pt->nnode();
6718 
6719  std::pair<Node*, Node*> tmp_pair = std::make_pair(
6720  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
6721 
6722  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
6723  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
6724 
6725  // Search for repeated nodes
6726  const unsigned n_done_nodes = done_nodes_pt.size();
6727  for (unsigned l = 0; l < n_done_nodes; l++)
6728  {
6729  if (tmp_pair == done_nodes_pt[l] ||
6730  tmp_pair_inverse == done_nodes_pt[l])
6731  {
6732  // Increase the number of repeated elements
6733  n_repeated_ele++;
6734  // Mark the element as repeated
6735  repeated = true;
6736  break;
6737  }
6738  }
6739 
6740  // Create new face element
6741  if (!repeated)
6742  {
6743  // Add the pair of nodes (edge) to the node dones
6744  done_nodes_pt.push_back(tmp_pair);
6745  // Add the element to the face elements
6746  face_el_pt.push_back(tmp_ele_pt);
6747  }
6748  else
6749  {
6750  // Free the repeated bulk element!!
6751  delete tmp_ele_pt;
6752  tmp_ele_pt = 0;
6753  }
6754 
6755  // Re-start
6756  repeated = false;
6757 
6758  } // for (e < nel)
6759  } // if (nel > 0)
6760 
6761  } // else (n_regions > 1)
6762 
6763  // Do not consider the repeated elements
6764  nele -= n_repeated_ele;
6765 
6766 #ifdef PARANOID
6767  if (nele != face_el_pt.size())
6768  {
6769  std::ostringstream error_message;
6770  error_message
6771  << "The independet counting of face elements (" << nele << ") for "
6772  << "boundary (" << b << ") is different\n"
6773  << "from the real number of face elements in the container ("
6774  << face_el_pt.size() << ")\n";
6775  //<< "Possible memory leak\n"
6776  throw OomphLibError(error_message.str(),
6777  "TriangleMesh::re_scale_re_assigned_initial_zeta_"
6778  "values_for_internal_boundary()",
6779  OOMPH_EXCEPTION_LOCATION);
6780  }
6781 #endif
6782 
6783  // ----------------------------------------------------------------
6784  // Second: Sort the face elements (to create segments), only
6785  // consider nonhalo elements
6786  // ----------------------------------------------------------------
6787 
6788  // Get the total number of nonhalo face elements
6789  const unsigned nnon_halo_face_elements = face_el_pt.size();
6790 
6791  // The vector of list to store the "segments" that compound the
6792  // boundary (segments may appear only in a distributed mesh)
6793  Vector<std::list<FiniteElement*>> segment_sorted_ele_pt;
6794 
6795  // Number of already sorted face elements
6796  unsigned nsorted_face_elements = 0;
6797 
6798  // Keep track of who's done
6799  std::map<FiniteElement*, bool> done_el;
6800 
6801  // Keep track of which element is inverted
6802  std::map<FiniteElement*, bool> is_inverted;
6803 
6804  // Iterate until all possible segments have been created
6805  while (nsorted_face_elements < nnon_halo_face_elements)
6806  {
6807  // The ordered list of face elements (in a distributed mesh a
6808  // collection of contiguous face elements define a segment)
6809  std::list<FiniteElement*> sorted_el_pt;
6810 
6811 #ifdef PARANOID
6812  // Select an initial element for the segment
6813  bool found_initial_face_element = false;
6814 #endif
6815 
6816  FiniteElement* ele_face_pt = 0;
6817 
6818  unsigned iface = 0;
6819  for (iface = 0; iface < nele; iface++)
6820  {
6821  ele_face_pt = face_el_pt[iface];
6822  // If not done then take it as initial face element
6823  if (!done_el[ele_face_pt])
6824  {
6825 #ifdef PARANOID
6826  found_initial_face_element = true;
6827 #endif
6828  nsorted_face_elements++;
6829  iface++; // The next element number
6830  sorted_el_pt.push_back(ele_face_pt);
6831  // Mark as done
6832  done_el[ele_face_pt] = true;
6833  break;
6834  }
6835  } // for (iface < nele)
6836 
6837 #ifdef PARANOID
6838  if (!found_initial_face_element)
6839  {
6840  std::ostringstream error_message;
6841  error_message
6842  << "Could not find an initial face element for the current segment\n";
6843  // << "----- Possible memory leak -----\n";
6844  throw OomphLibError(error_message.str(),
6845  "TriangleMesh::re_scale_re_assigned_initial_zeta_"
6846  "values_for_internal_boundary()",
6847  OOMPH_EXCEPTION_LOCATION);
6848  }
6849 #endif
6850 
6851  // Number of nodes
6852  const unsigned nnod = ele_face_pt->nnode();
6853 
6854  // Left and rightmost nodes (the left and right nodes of the
6855  // current face element)
6856  Node* left_node_pt = ele_face_pt->node_pt(0);
6857  Node* right_node_pt = ele_face_pt->node_pt(nnod - 1);
6858 
6859  // Continue iterating if a new face element has been added to the
6860  // list
6861  bool face_element_added = false;
6862 
6863  // While a new face element has been added to the set of sorted
6864  // face elements then re-iterate
6865  do
6866  {
6867  // Start from the next face element since we have already added
6868  // the previous one as the initial face element (any previous
6869  // face element had to be added on previous iterations)
6870  for (unsigned iiface = iface; iiface < nele; iiface++)
6871  {
6872  // Re-start flag
6873  face_element_added = false;
6874 
6875  // Get the candidate element
6876  ele_face_pt = face_el_pt[iiface];
6877 
6878  // Check that the candidate element has not been done
6879  if (!(done_el[ele_face_pt]))
6880  {
6881  // Get the left and right nodes of the current element
6882  Node* local_left_node_pt = ele_face_pt->node_pt(0);
6883  Node* local_right_node_pt = ele_face_pt->node_pt(nnod - 1);
6884 
6885  // New element fits at the left of segment and is not inverted
6886  if (left_node_pt == local_right_node_pt)
6887  {
6888  left_node_pt = local_left_node_pt;
6889  sorted_el_pt.push_front(ele_face_pt);
6890  is_inverted[ele_face_pt] = false;
6891  face_element_added = true;
6892  }
6893  // New element fits at the left of segment and is inverted
6894  else if (left_node_pt == local_left_node_pt)
6895  {