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