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