multi_domain.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// Templated multi-domain functions
27
28// Include guards to prevent multiple inclusion of the header
29#ifndef OOMPH_MULTI_DOMAIN_CC
30#define OOMPH_MULTI_DOMAIN_CC
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37// Oomph-lib headers
38#include "geom_objects.h"
39#include "problem.h"
40#include "shape.h"
41
42#include "mesh.h"
44#include "algebraic_elements.h"
46#include "Qelements.h"
48#include "multi_domain.h"
50
51// Needed to check if elements have nonuniformly spaced nodes
52#include "refineable_elements.h"
53#include "Qspectral_elements.h"
54
55namespace oomph
56{
57 /// / Templated helper functions for multi-domain methods using locate_zeta
58
59
60 //============================================================================
61 /// Identify the \c FaceElements (stored in the mesh pointed to by
62 /// \c face_mesh_pt) that are adjacent to the bulk elements next to the
63 /// \c boundary_in_bulk_mesh -th boundary of the mesh pointed to by
64 /// \c bulk_mesh_pt. The \c FaceElements must be derived
65 /// from the \c ElementWithExternalElement base class and the adjacent
66 /// bulk elements are stored as their external elements.
67 ///
68 /// This is the vector-based version which deals with multiple bulk
69 /// mesh boundaries at the same time.
70 //============================================================================
71 template<class BULK_ELEMENT, unsigned DIM>
73 Problem* problem_pt,
74 Vector<unsigned>& boundary_in_bulk_mesh,
75 Mesh* const& bulk_mesh_pt,
76 Vector<Mesh*>& face_mesh_pt,
77 const unsigned& interaction)
78 {
79 unsigned n_mesh = boundary_in_bulk_mesh.size();
80
81#ifdef PARANOID
82 // Check sizes match
83 if (boundary_in_bulk_mesh.size() != face_mesh_pt.size())
84 {
85 std::ostringstream error_message;
86 error_message << "Sizes of vector of boundary ids in bulk mesh ("
87 << boundary_in_bulk_mesh.size()
88 << ") and vector of pointers\n"
89 << "to FaceElements (" << face_mesh_pt.size()
90 << " doesn't match.\n";
91 throw OomphLibError(
92 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
93 }
94#endif
95
96 // Create face meshes adjacent to the bulk mesh's b-th boundary.
97 // Each face mesh consists of FaceElements that may also be
98 // interpreted as GeomObjects
99 Vector<Mesh*> bulk_face_mesh_pt(n_mesh);
100
101 // Loop over all meshes
102 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
103 {
104 bulk_face_mesh_pt[i_mesh] = new Mesh;
105 bulk_mesh_pt
106 ->template build_face_mesh<BULK_ELEMENT, FaceElementAsGeomObject>(
107 boundary_in_bulk_mesh[i_mesh], bulk_face_mesh_pt[i_mesh]);
108
109 // Loop over these new face elements and tell them the boundary number
110 // from the bulk mesh -- this is required to they can
111 // get access to the boundary coordinates!
112 unsigned n_face_element = bulk_face_mesh_pt[i_mesh]->nelement();
113 for (unsigned e = 0; e < n_face_element; e++)
114 {
115 // Cast the element pointer to the correct thing!
118 bulk_face_mesh_pt[i_mesh]->element_pt(e));
119
120 // Set bulk boundary number
121 el_pt->set_boundary_number_in_bulk_mesh(boundary_in_bulk_mesh[i_mesh]);
122
123 // Doc?
124 if (Doc_boundary_coordinate_file.is_open())
125 {
126 Vector<double> s(DIM - 1);
127 Vector<double> zeta(DIM - 1);
128 Vector<double> x(DIM);
129 unsigned n_plot = 5;
131
132 // Loop over plot points
133 unsigned num_plot_points = el_pt->nplot_points(n_plot);
134 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
135 {
136 // Get local coordinates of plot point
137 el_pt->get_s_plot(iplot, n_plot, s);
138 el_pt->interpolated_zeta(s, zeta);
139 el_pt->interpolated_x(s, x);
140 for (unsigned i = 0; i < DIM; i++)
141 {
142 Doc_boundary_coordinate_file << x[i] << " ";
143 }
144 for (unsigned i = 0; i < DIM - 1; i++)
145 {
146 Doc_boundary_coordinate_file << zeta[i] << " ";
147 }
148 Doc_boundary_coordinate_file << std::endl;
149 }
151 n_plot);
152 }
153 }
154
155 // Now sort the elements based on the boundary coordinates.
156 // This may allow a faster implementation of the locate_zeta
157 // function for the MeshAsGeomObject representation of this
158 // mesh, but also creates a unique ordering of the elements
159 std::sort(bulk_face_mesh_pt[i_mesh]->element_pt().begin(),
160 bulk_face_mesh_pt[i_mesh]->element_pt().end(),
162 } // end of loop over meshes
163
164
165 // Setup the interactions for this problem using the surface mesh
166 // on the fluid mesh. The interaction parameter in this instance is
167 // given by the "face" parameter.
169 BULK_ELEMENT,
171 problem_pt, face_mesh_pt, bulk_mesh_pt, bulk_face_mesh_pt, interaction);
172
173
174 // Loop over all meshes to clean up
175 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
176 {
177 unsigned n_face_element = bulk_face_mesh_pt[i_mesh]->nelement();
178
179 // The MeshAsGeomObject has already been deleted (in set_external_storage)
180
181 // Must be careful with the FaceMesh, because we cannot delete the nodes
182 // Loop over the elements backwards (paranoid!) and delete them
183 for (unsigned e = n_face_element; e > 0; e--)
184 {
185 delete bulk_face_mesh_pt[i_mesh]->element_pt(e - 1);
186 bulk_face_mesh_pt[i_mesh]->element_pt(e - 1) = 0;
187 }
188 // Now clear all element and node storage (should maybe fine-grain this)
189 bulk_face_mesh_pt[i_mesh]->flush_element_and_node_storage();
190
191 // Finally delete the mesh
192 delete bulk_face_mesh_pt[i_mesh];
193
194 } // end of loop over meshes
195 }
196
197
198 //========================================================================
199 /// Identify the \c FaceElements (stored in the mesh pointed to by
200 /// \c face_mesh_pt) that are adjacent to the bulk elements next to the
201 /// \c boundary_in_bulk_mesh -th boundary of the mesh pointed to by
202 /// \c bulk_mesh_pt. The \c FaceElements must be derived
203 /// from the \c ElementWithExternalElement base class and the adjacent
204 /// bulk elements are stored as their external elements.
205 //========================================================================
206 template<class BULK_ELEMENT, unsigned DIM>
208 Problem* problem_pt,
209 const unsigned& boundary_in_bulk_mesh,
210 Mesh* const& bulk_mesh_pt,
211 Mesh* const& face_mesh_pt,
212 const unsigned& interaction)
213 {
214 // Convert to vector-based storage
215 Vector<unsigned> boundary_in_bulk_mesh_vect(1);
216 boundary_in_bulk_mesh_vect[0] = boundary_in_bulk_mesh;
217 Vector<Mesh*> face_mesh_pt_vect(1);
218 face_mesh_pt_vect[0] = face_mesh_pt;
219
220 // Call vector-based version
221 setup_bulk_elements_adjacent_to_face_mesh<BULK_ELEMENT, DIM>(
222 problem_pt,
223 boundary_in_bulk_mesh_vect,
224 bulk_mesh_pt,
225 face_mesh_pt_vect,
226 interaction);
227 }
228
229
230 //========================================================================
231 /// Set up the two-way multi-domain interactions for the
232 /// problem pointed to by \c problem_pt.
233 /// Use this for cases where first_mesh_pt and second_mesh_pt
234 /// occupy the same physical space and are populated by
235 /// ELEMENT_0 and ELEMENT_1 respectively, and are combined to solve
236 /// a single problem. The elements in two meshes interact both ways
237 /// the elements in each mesh act as "external elements" for the
238 /// elements in the "other" mesh. The interaction indices allow the
239 /// specification of which interaction we're setting up in the two
240 /// meshes. They default to zero, which is appropriate if there's
241 /// only a single interaction.
242 //========================================================================
243 template<class ELEMENT_0, class ELEMENT_1>
245 Problem* problem_pt,
246 Mesh* const& first_mesh_pt,
247 Mesh* const& second_mesh_pt,
248 const unsigned& first_interaction,
249 const unsigned& second_interaction)
250 {
251 // Delete all the current external halo(ed) element and node storage
252 first_mesh_pt->delete_all_external_storage();
253 second_mesh_pt->delete_all_external_storage();
254
255 // Call setup_multi_domain_interaction in both directions
256 setup_multi_domain_interaction<ELEMENT_1>(
257 problem_pt, first_mesh_pt, second_mesh_pt, first_interaction);
258
259 setup_multi_domain_interaction<ELEMENT_0>(
260 problem_pt, second_mesh_pt, first_mesh_pt, second_interaction);
261 }
262
263
264 //========================================================================
265 /// Function to set up the one-way multi-domain interaction for
266 /// problems where the meshes pointed to by \c mesh_pt and \c external_mesh_pt
267 /// occupy the same physical space, and the elements in \c external_mesh_pt
268 /// act as "external elements" for the \c ElementWithExternalElements
269 /// in \c mesh_pt (but not vice versa):
270 /// - \c mesh_pt points to the mesh of ElemenWithExternalElements for which
271 /// the interaction is set up.
272 /// - \c external_mesh_pt points to the mesh that contains the elements
273 /// of type EXT_ELEMENT that act as "external elements" for the
274 /// \c ElementWithExternalElements in \c mesh_pt.
275 /// - The interaction_index parameter defaults to zero and must be otherwise
276 /// set by the user if there is more than one mesh that provides sources
277 /// for the Mesh pointed to by mesh_pt.
278 //========================================================================
279 template<class EXT_ELEMENT>
281 Problem* problem_pt,
282 Mesh* const& mesh_pt,
283 Mesh* const& external_mesh_pt,
284 const unsigned& interaction_index)
285 {
286 // Bulk elements must not be external elements in this case
288
289 // Call the auxiliary function with GEOM_OBJECT=EXT_ELEMENT
290 // and EL_DIM_EUL=EL_DIM_LAG=dimension returned from helper function
291 get_dim_helper(problem_pt, mesh_pt, external_mesh_pt);
292
293 if (Dim > 3)
294 {
295 std::ostringstream error_stream;
296 error_stream << "The elements within the two interacting meshes have a\n"
297 << "dimension not equal to 1, 2 or 3.\n"
298 << "The multi-domain method will not work in this case.\n"
299 << "The dimension is: " << Dim << "\n";
300 throw OomphLibError(
301 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
302 }
303
304 // Wrapper for each dimension (template parameter)
305 aux_setup_multi_domain_interaction<EXT_ELEMENT, EXT_ELEMENT>(
306 problem_pt, mesh_pt, external_mesh_pt, interaction_index);
307 }
308
309
310 //========================================================================
311 /// Function to set up the one-way multi-domain interaction for
312 /// FSI-like problems.
313 /// - \c mesh_pt points to the mesh of \c ElemenWithExternalElements for which
314 /// the interaction is set up. In an FSI example, this mesh would contain
315 /// the \c FSIWallElements (either beam/shell elements or the
316 /// \c FSISolidTractionElements that apply the traction to
317 /// a "bulk" solid mesh that is loaded by the fluid.)
318 /// - \c external_mesh_pt points to the mesh that contains the elements
319 /// of type EXT_ELEMENT that provide the "source" for the
320 /// \c ElementWithExternalElements. In an FSI example, this
321 /// mesh would contain the "bulk" fluid elements.
322 /// - \c external_face_mesh_pt points to the mesh of \c FaceElements
323 /// attached to the \c external_mesh_pt. The mesh pointed to by
324 /// \c external_face_mesh_pt has the same dimension as \c mesh_pt.
325 /// The elements contained in \c external_face_mesh_pt are of type
326 /// FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements
327 /// are usually the \c FaceElementAsGeomObjects (templated by the
328 /// type of the "bulk" fluid elements to which they are attached)
329 /// that define the FSI boundary of the fluid domain.
330 /// - The interaction_index parameter defaults to zero and must otherwise be
331 /// set by the user if there is more than one mesh that provides "external
332 /// elements" for the Mesh pointed to by mesh_pt (e.g. in the case
333 /// when a beam or shell structure is loaded by fluid from both sides.)
334 //========================================================================
335 template<class EXT_ELEMENT, class FACE_ELEMENT_GEOM_OBJECT>
337 Problem* problem_pt,
338 Mesh* const& mesh_pt,
339 Mesh* const& external_mesh_pt,
340 Mesh* const& external_face_mesh_pt,
341 const unsigned& interaction_index)
342 {
343 // Bulk elements must be external elements in this case
345
346 // Call the auxiliary routine with GEOM_OBJECT=FACE_ELEMENT_GEOM_OBJECT
347 // and EL_DIM_LAG=Dim, EL_DIM_EUL=Dim+1
348 get_dim_helper(problem_pt, mesh_pt, external_face_mesh_pt);
349
350 if (Dim > 2)
351 {
352 std::ostringstream error_stream;
353 error_stream << "The elements within the two interacting meshes have a\n"
354 << "dimension not equal to 1 or 2.\n"
355 << "The multi-domain method will not work in this case.\n"
356 << "The dimension is: " << Dim << "\n";
357 throw OomphLibError(
358 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
359 }
360
361 // Call the function that actually does the work
362 aux_setup_multi_domain_interaction<EXT_ELEMENT, FACE_ELEMENT_GEOM_OBJECT>(
363 problem_pt,
364 mesh_pt,
365 external_mesh_pt,
366 interaction_index,
367 external_face_mesh_pt);
368 }
369
370
371 //========================================================================
372 /// Function to set up the one-way multi-domain interaction for
373 /// FSI-like problems.
374 /// - \c mesh_pt points to the mesh of \c ElemenWithExternalElements for which
375 /// the interaction is set up. In an FSI example, this mesh would contain
376 /// the \c FSIWallElements (either beam/shell elements or the
377 /// \c FSISolidTractionElements that apply the traction to
378 /// a "bulk" solid mesh that is loaded by the fluid.)
379 /// - \c external_mesh_pt points to the mesh that contains the elements
380 /// of type EXT_ELEMENT that provide the "source" for the
381 /// \c ElementWithExternalElements. In an FSI example, this
382 /// mesh would contain the "bulk" fluid elements.
383 /// - \c external_face_mesh_pt points to the mesh of \c FaceElements
384 /// attached to the \c external_mesh_pt. The mesh pointed to by
385 /// \c external_face_mesh_pt has the same dimension as \c mesh_pt.
386 /// The elements contained in \c external_face_mesh_pt are of type
387 /// FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements
388 /// are usually the \c FaceElementAsGeomObjects (templated by the
389 /// type of the "bulk" fluid elements to which they are attached)
390 /// that define the FSI boundary of the fluid domain.
391 /// - The interaction_index parameter defaults to zero and must otherwise be
392 /// set by the user if there is more than one mesh that provides "external
393 /// elements" for the Mesh pointed to by mesh_pt (e.g. in the case
394 /// when a beam or shell structure is loaded by fluid from both sides.)
395 /// .
396 /// Vector-based version operates simultaneously on the meshes contained
397 /// in the vectors.
398 //========================================================================
399 template<class EXT_ELEMENT, class FACE_ELEMENT_GEOM_OBJECT>
401 Problem* problem_pt,
402 const Vector<Mesh*>& mesh_pt,
403 Mesh* const& external_mesh_pt,
404 const Vector<Mesh*>& external_face_mesh_pt,
405 const unsigned& interaction_index)
406 {
407 // How many meshes do we have?
408 unsigned n_mesh = mesh_pt.size();
409
410#ifdef PARANOID
411 if (external_face_mesh_pt.size() != n_mesh)
412 {
413 std::ostringstream error_stream;
414 error_stream << "Sizes of external_face_mesh_pt [ "
415 << external_face_mesh_pt.size() << " ] and "
416 << "mesh_pt [ " << n_mesh << " ] don't match.\n";
417 throw OomphLibError(
418 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
419 }
420#endif
421
422 // Bail out?
423 if (n_mesh == 0) return;
424
425 // Bulk elements must be external elements in this case
427
428 // Call the auxiliary routine with GEOM_OBJECT=FACE_ELEMENT_GEOM_OBJECT
429 // and EL_DIM_LAG=Dim, EL_DIM_EUL=Dim+1. Use first mesh only.
430 get_dim_helper(problem_pt, mesh_pt[0], external_face_mesh_pt[0]);
431
432
433#ifdef PARANOID
434 // Check consitency
435 unsigned old_dim = Dim;
436 for (unsigned i = 1; i < n_mesh; i++)
437 {
438 // Set up Dim again
439 get_dim_helper(problem_pt, mesh_pt[i], external_face_mesh_pt[i]);
440
441 if (Dim != old_dim)
442 {
443 std::ostringstream error_stream;
444 error_stream << "Inconsistency: Mesh " << i << " has Dim=" << Dim
445 << "while mesh 0 has Dim=" << old_dim << std::endl;
446 throw OomphLibError(
447 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
448 }
449 }
450#endif
451
452 if (Dim > 2)
453 {
454 std::ostringstream error_stream;
455 error_stream << "The elements within the two interacting meshes have a\n"
456 << "dimension not equal to 1 or 2.\n"
457 << "The multi-domain method will not work in this case.\n"
458 << "The dimension is: " << Dim << "\n";
459 throw OomphLibError(
460 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
461 }
462
463 // Now do the actual work for all meshes simultaneously
464 aux_setup_multi_domain_interaction<EXT_ELEMENT, FACE_ELEMENT_GEOM_OBJECT>(
465 problem_pt,
466 mesh_pt,
467 external_mesh_pt,
468 interaction_index,
469 external_face_mesh_pt);
470 }
471
472
473 //========================================================================
474 /// This routine calls the locate_zeta routine (simultaneously on each
475 /// processor for each individual processor's element set if necessary)
476 /// and sets up the external (halo) element and node storage as
477 /// necessary. The locate_zeta procedure here works for all multi-domain
478 /// problems where either two meshes occupy the same physical space but have
479 /// differing element types (e.g. a Boussinesq convection problem where
480 /// AdvectionDiffusion elements interact with Navier-Stokes type elements)
481 /// or two meshes interact along some boundary of the external mesh,
482 /// represented by a "face mesh", such as an FSI problem.
483 //========================================================================
484 template<class EXT_ELEMENT, class GEOM_OBJECT>
486 Problem* problem_pt,
487 Mesh* const& mesh_pt,
488 Mesh* const& external_mesh_pt,
489 const unsigned& interaction_index,
490 Mesh* const& external_face_mesh_pt)
491 {
492 // Convert to vector-based storage
493 Vector<Mesh*> mesh_pt_vector(1);
494 mesh_pt_vector[0] = mesh_pt;
495 Vector<Mesh*> external_face_mesh_pt_vector(1);
496 external_face_mesh_pt_vector[0] = external_face_mesh_pt;
497
498 // Call vector-based version
499 aux_setup_multi_domain_interaction<EXT_ELEMENT, GEOM_OBJECT>(
500 problem_pt,
501 mesh_pt_vector,
502 external_mesh_pt,
503 interaction_index,
504 external_face_mesh_pt_vector);
505
506 } // end of aux_setup_multi_domain_interaction
507
508
509 //========================================================================
510 /// This routine calls the locate_zeta routine (simultaneously on each
511 /// processor for each individual processor's element set if necessary)
512 /// and sets up the external (halo) element and node storage as
513 /// necessary. The locate_zeta procedure here works for all multi-domain
514 /// problems where either two meshes occupy the same physical space but have
515 /// differing element types (e.g. a Boussinesq convection problem where
516 /// AdvectionDiffusion elements interact with Navier-Stokes type elements)
517 /// or two meshes interact along some boundary of the external mesh,
518 /// represented by a "face mesh", such as an FSI problem.
519 ///
520 /// Vector-based version operates simultaneously on the meshes contained
521 /// in the vectors.
522 //========================================================================
523 template<class EXT_ELEMENT, class GEOM_OBJECT>
525 Problem* problem_pt,
526 const Vector<Mesh*>& mesh_pt,
527 Mesh* const& external_mesh_pt,
528 const unsigned& interaction_index,
529 const Vector<Mesh*>& external_face_mesh_pt)
530 {
531 // How many meshes do we have?
532 unsigned n_mesh = mesh_pt.size();
533
534#ifdef PARANOID
535 if (external_face_mesh_pt.size() != n_mesh)
536 {
537 std::ostringstream error_stream;
538 error_stream << "Sizes of external_face_mesh_pt [ "
539 << external_face_mesh_pt.size() << " ] and "
540 << "mesh_pt [ " << n_mesh << " ] don't match.\n";
541 throw OomphLibError(
542 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
543 }
544#endif
545
546 // Bail out?
547 if (n_mesh == 0) return;
548
549 // Multi-domain setup will not work for elements with
550 // nonuniformly spaced nodes
551 // Must check type of elements in the mesh and in the external mesh
552 //(assume element type is the same for all elements in each mesh)
553
554#ifdef PARANOID
555
556 // Pointer to first element in external mesh
557 GeneralisedElement* ext_el_pt_0 = 0;
558 if (external_mesh_pt->nelement() != 0)
559 {
560 ext_el_pt_0 = external_mesh_pt->element_pt(0);
561 }
562
563 // Loop over all meshes
564 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
565 {
566 // Get pointer to first element in each mesh
567 GeneralisedElement* el_pt_0 = 0;
568 if (mesh_pt[i_mesh]->nelement() != 0)
569 {
570 el_pt_0 = mesh_pt[i_mesh]->element_pt(0);
571 }
572
573 // Check they are not spectral elements
574 if (dynamic_cast<SpectralElement*>(el_pt_0) != 0 ||
575 dynamic_cast<SpectralElement*>(ext_el_pt_0) != 0)
576 {
577 throw OomphLibError(
578 "Multi-domain setup does not work with spectral elements.",
579 OOMPH_CURRENT_FUNCTION,
580 OOMPH_EXCEPTION_LOCATION);
581 }
582
583 // Check they are not hp-refineable elements
584 if (dynamic_cast<PRefineableElement*>(el_pt_0) != 0 ||
585 dynamic_cast<PRefineableElement*>(ext_el_pt_0) != 0)
586 {
587 throw OomphLibError(
588 "Multi-domain setup does not work with hp-refineable elements.",
589 OOMPH_CURRENT_FUNCTION,
590 OOMPH_EXCEPTION_LOCATION);
591 }
592 } // end over initial loop over meshes
593
594#endif
595
596
597#ifdef OOMPH_HAS_MPI
598 // Storage for number of processors and my rank
599 int n_proc = problem_pt->communicator_pt()->nproc();
600 int my_rank = problem_pt->communicator_pt()->my_rank();
601#endif
602
603 // Timing
604 double t_start = 0.0;
605 double t_end = 0.0;
606 double t_local = 0.0;
607 double t_set = 0.0;
608 double t_locate = 0.0;
609 double t_spiral_start = 0.0;
610#ifdef OOMPH_HAS_MPI
611 double t_loop_start = 0.0;
612 double t_sendrecv = 0.0;
613 double t_missing = 0.0;
614 double t_send_info = 0.0;
615 double t_create_halo = 0.0;
616#endif
617
618 if (Doc_timings)
619 {
620 t_start = TimingHelpers::timer();
621 }
622
623 // Initialize number of zeta coordinates not found yet
624 unsigned n_zeta_not_found = 0;
625
626 // Geometric objects used to represent the external (face) meshes
627 Vector<MeshAsGeomObject*> mesh_geom_obj_pt(n_mesh, 0);
628
629#ifdef PARANOID
630
631 // Initialise lagrangian dimension of element (test only)
632 unsigned el_dim_lag = 0;
633
634#endif
635
636 // Create mesh as geom objects for all meshes
637 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
638 {
639 // Are bulk elements used as external elements?
641 {
642 // Fix this when required
643 if (n_mesh != 1)
644 {
645 std::ostringstream error_stream;
646 error_stream
647 << "Sorry I currently can't deal with non-bulk external elements\n"
648 << "in multi-domain setup for multiple meshes.\n"
649 << "The functionality should be easy to implement now that you\n"
650 << "have a test case. If you're not willinig to do this, call\n"
651 << "the multi-domain setup mesh-by-mesh instead (though this can\n"
652 << "be costly in parallel because of global comms. \n";
653 throw OomphLibError(error_stream.str(),
654 OOMPH_CURRENT_FUNCTION,
655 OOMPH_EXCEPTION_LOCATION);
656 }
657
658 // Set the geometric object from the external mesh
659 mesh_geom_obj_pt[0] = new MeshAsGeomObject(external_mesh_pt);
660 }
661 else
662 {
663 // Set the geometric object from the external face mesh argument
664 mesh_geom_obj_pt[i_mesh] =
665 new MeshAsGeomObject(external_face_mesh_pt[i_mesh]);
666 }
667
668#ifdef PARANOID
669 unsigned old_el_dim_lag = el_dim_lag;
670
671 // Set lagrangian dimension of element
672 el_dim_lag = mesh_geom_obj_pt[i_mesh]->nlagrangian();
673
674 // Check consistency
675 if (i_mesh > 0)
676 {
677 if (el_dim_lag != old_el_dim_lag)
678 {
679 std::ostringstream error_stream;
680 error_stream << "Lagrangian dimensions of elements don't match \n "
681 << "between meshes: " << el_dim_lag << " "
682 << old_el_dim_lag << "\n";
683 throw OomphLibError(error_stream.str(),
684 OOMPH_CURRENT_FUNCTION,
685 OOMPH_EXCEPTION_LOCATION);
686 }
687 }
688#endif
689
690
691 } // end of loop over meshes
692
693 double t_setup_lookups = 0.0;
694 if (Doc_timings)
695 {
696 t_set = TimingHelpers::timer();
697 oomph_info << "CPU for creation of MeshAsGeomObjects and bin structure: "
698 << t_set - t_start << std::endl;
699 t_setup_lookups = TimingHelpers::timer();
700 }
701
702 // Total number of integration points
703 unsigned tot_int = 0;
704
705 // Counter for total number of elements (in flat packed order)
706 unsigned e_count = 0;
707
708 // Loop over all meshes to get total number of elements
709 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
710 {
711 e_count += mesh_pt[i_mesh]->nelement();
712 }
713 External_element_located.resize(e_count);
714
715 // Reset counter for elements in flat packed storage
716 e_count = 0;
717
718 // Loop over all meshes
719 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
720 {
721 // Loop over (this processor's) elements and set lookup array
722 unsigned n_element = mesh_pt[i_mesh]->nelement();
723 for (unsigned e = 0; e < n_element; e++)
724 {
725 // Zero-sized vector means its a halo
726 External_element_located[e_count].resize(0);
728 dynamic_cast<ElementWithExternalElement*>(
729 mesh_pt[i_mesh]->element_pt(e));
730
731#ifdef OOMPH_HAS_MPI
732 // We're not setting up external elements for halo elements
733 if (!el_pt->is_halo())
734#endif
735 {
736 // We need to allocate storage for the external elements
737 // within the element. Memory will actually only be
738 // allocated the first time this function is called for
739 // each element, or if the number of interactions or integration
740 // points within the element has changed.
742
743 // Clear any previous allocation
744 unsigned n_intpt = el_pt->integral_pt()->nweight();
745 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
746 {
747 el_pt->external_element_pt(interaction_index, ipt) = 0;
748 }
749
750 External_element_located[e_count].resize(n_intpt);
751 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
752 {
753 External_element_located[e_count][ipt] = 0;
754 tot_int++;
755 }
756 }
757 // next element
758 e_count++;
759 }
760 } // end of loop over meshes
761
762 if (Doc_timings)
763 {
764 double t = TimingHelpers::timer();
766 << "CPU for setup of lookup schemes for located elements/coords: "
767 << t - t_setup_lookups << std::endl;
768 }
769
770 // Initialise maximum spiral level within the cartesian bin structure
771 // Used to terminate spiraling for non-refineable bin
772 unsigned n_max_level = 0;
773
774#ifdef OOMPH_HAS_MPI
775 unsigned max_level_reached = 1;
776#endif
777
778 // Max. number of sample points -- used to decide on termination of
779 // "spiraling"
780 unsigned max_n_sample_points_of_sample_point_containers = 0;
781
782 // Loop over all meshes
783 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
784 {
785 if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
787 {
788 RefineableBinArray* bin_array_pt = dynamic_cast<RefineableBinArray*>(
789 mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
790
791 bin_array_pt
793 bin_array_pt
795 bin_array_pt
797
798 unsigned nsp =
800 if (nsp > max_n_sample_points_of_sample_point_containers)
801 {
802 max_n_sample_points_of_sample_point_containers = nsp;
803 }
804
805
806#ifdef OOMPH_HAS_MPI
807 // If the mesh has been distributed we want the max. number
808 // of sample points across all processors
809 if (problem_pt->communicator_pt()->nproc() > 1)
810 {
811 unsigned local_max_n_sample_points_of_sample_point_containers =
812 max_n_sample_points_of_sample_point_containers;
813
814 // Get maximum over all processors
815 MPI_Allreduce(&local_max_n_sample_points_of_sample_point_containers,
816 &max_n_sample_points_of_sample_point_containers,
817 1,
818 MPI_UNSIGNED,
819 MPI_MAX,
820 problem_pt->communicator_pt()->mpi_comm());
821 }
822#endif
823 }
824 else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
826 {
827 NonRefineableBinArray* bin_array_pt =
828 dynamic_cast<NonRefineableBinArray*>(
829 mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
830
831 // Initialise spiral levels
832 bin_array_pt->current_min_spiral_level() = 0;
833 bin_array_pt->current_max_spiral_level() =
834 bin_array_pt->n_spiral_chunk() - 1;
835
836 // Find maximum spiral level within the cartesian bin structure
837 n_max_level = bin_array_pt->max_bin_dimension();
838
839 // Limit it
840 if (bin_array_pt->current_max_spiral_level() > n_max_level)
841 {
842 bin_array_pt->current_max_spiral_level() = n_max_level - 1;
843 }
844 }
845#ifdef OOMPH_HAS_CGAL
846 // CGAL
847 else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
849 {
850 CGALSamplePointContainer* bin_array_pt =
851 dynamic_cast<CGALSamplePointContainer*>(
852 mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
853 bin_array_pt
855 bin_array_pt
857 bin_array_pt
859
860 unsigned nsp =
862 if (nsp > max_n_sample_points_of_sample_point_containers)
863 {
864 max_n_sample_points_of_sample_point_containers = nsp;
865 }
866
867
868#ifdef OOMPH_HAS_MPI
869 // If the mesh has been distributed we want the max. number
870 // of sample points across all processors
871 if (problem_pt->communicator_pt()->nproc() > 1)
872 {
873 unsigned local_max_n_sample_points_of_sample_point_containers =
874 max_n_sample_points_of_sample_point_containers;
875
876 // Get maximum over all processors
877 MPI_Allreduce(&local_max_n_sample_points_of_sample_point_containers,
878 &max_n_sample_points_of_sample_point_containers,
879 1,
880 MPI_UNSIGNED,
881 MPI_MAX,
882 problem_pt->communicator_pt()->mpi_comm());
883 }
884#endif
885 }
886#endif // cgal
887 }
888
889
890 // Storage for info about coordinate location
891 Vector<double> percentage_coords_located_locally;
892 Vector<double> percentage_coords_located_elsewhere;
893
894 // Loop over "spirals/levels" away from the current position
895 // Note: All meshes go through their spirals simultaneously;
896 // read out spiral level from first one
897 unsigned i_level = 0;
898 bool has_not_reached_max_level_of_search = true;
899 while (has_not_reached_max_level_of_search)
900 {
901 // Record time at start of spiral loop
902 if (Doc_timings)
903 {
904 t_spiral_start = TimingHelpers::timer();
905 }
906
907 // Perform locate_zeta locally first! This looks locally for
908 // all not-yet-located zetas within the current spiral range.
910 mesh_pt, external_mesh_pt, mesh_geom_obj_pt, interaction_index);
911
912 // Store stats about successful locates for reporting later
913 if (Doc_stats)
914 {
915 unsigned count_locates = 0;
916 unsigned n_ext_loc = External_element_located.size();
917 for (unsigned e = 0; e < n_ext_loc; e++)
918 {
919 unsigned n_intpt = External_element_located[e].size();
920 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
921 {
922 count_locates += External_element_located[e][ipt];
923 }
924 }
925
926 // Store percentage of integration points successfully located.
927 // Only assign if we had anything to allocte, otherwise 100%
928 // (default assignment; see above) is correct
929 if (tot_int != 0)
930 {
931 percentage_coords_located_locally.push_back(
932 100.0 * double(count_locates) / double(tot_int));
933 }
934 else
935 {
936 // Had none to find so we found them all!
937 percentage_coords_located_locally.push_back(100.0);
938 }
939 }
940
941
942 // Now test whether anything needs to be broadcast elsewhere
943 // (i.e. were there any failures in the locate method above?)
944 // If there are, then the zetas for these failures need to be
945 // broadcast...
946
947 // How many zetas have we failed to find? [Note: Array is padded
948 // by Dim padded entries (DBL_MAX) for each mesh]
949 n_zeta_not_found =
951
952 if (Doc_timings)
953 {
954 t_local = TimingHelpers::timer();
955 oomph_info << "CPU for local location of zeta coordinate [spiral level "
956 << i_level << "]: " << t_local - t_spiral_start << std::endl
957 << "Number of missing zetas: " << n_zeta_not_found
958 << std::endl;
959 }
960
961
962#ifdef OOMPH_HAS_MPI
963 // Only perform the reduction operation if there's more than one process
964 if (problem_pt->communicator_pt()->nproc() > 1)
965 {
966 unsigned count_local_zetas = n_zeta_not_found;
967 MPI_Allreduce(&count_local_zetas,
968 &n_zeta_not_found,
969 1,
970 MPI_UNSIGNED,
971 MPI_SUM,
972 problem_pt->communicator_pt()->mpi_comm());
973 }
974
975 // If we have missing zetas on any process
976 // and the problem is distributed, we need to locate elsewhere
977 if ((n_zeta_not_found != 0) &&
978 (problem_pt->problem_has_been_distributed()))
979 {
980 // Timings
981 double t_sendrecv_min = DBL_MAX;
982 double t_sendrecv_max = -DBL_MAX;
983 double t_sendrecv_tot = 0.0;
984
985 double t_missing_min = DBL_MAX;
986 double t_missing_max = -DBL_MAX;
987 double t_missing_tot = 0.0;
988
989 double t_send_info_min = DBL_MAX;
990 double t_send_info_max = -DBL_MAX;
991 double t_send_info_tot = 0.0;
992
993 double t_create_halo_min = DBL_MAX;
994 double t_create_halo_max = -DBL_MAX;
995 double t_create_halo_tot = 0.0;
996
997 // Start ring communication: Loop (number of processes - 1)
998 // starting from 1. The variable iproc represents the "distance" from
999 // the current process to the process for which it is attempting
1000 // to locate an element for the current set of not-yet-located
1001 // zeta coordinates
1002 unsigned ring_count = 0;
1003 for (int iproc = 1; iproc < n_proc; iproc++)
1004 {
1005 // Record time at start of loop
1006 if (Doc_timings)
1007 {
1008 t_loop_start = TimingHelpers::timer();
1009 }
1010
1011 // Send the zeta values you haven't found to the
1012 // next process, receive from the previous process:
1013 // (Padded) Flat_packed_zetas_not_found_locally are sent
1014 // to next processor where they are received as
1015 // (padded) Received_flat_packed_zetas_to_be_found.
1017
1018 if (Doc_timings)
1019 {
1020 ring_count++;
1021 t_sendrecv = TimingHelpers::timer();
1022 t_sendrecv_max =
1023 std::max(t_sendrecv_max, t_sendrecv - t_loop_start);
1024 t_sendrecv_min =
1025 std::min(t_sendrecv_min, t_sendrecv - t_loop_start);
1026 t_sendrecv_tot += (t_sendrecv - t_loop_start);
1027 }
1028
1029 // Perform the locate_zeta for the new set of zetas on this process
1031 iproc, external_mesh_pt, problem_pt, mesh_geom_obj_pt);
1032
1033 if (Doc_timings)
1034 {
1035 t_missing = TimingHelpers::timer();
1036 t_missing_max = std::max(t_missing_max, t_missing - t_sendrecv);
1037 t_missing_min = std::min(t_missing_min, t_missing - t_sendrecv);
1038 t_missing_tot += (t_missing - t_sendrecv);
1039 }
1040
1041 // Send any located coordinates back to the correct process, and
1042 // prepare to send on to the next process if necessary
1043 send_and_receive_located_info(iproc, external_mesh_pt, problem_pt);
1044
1045 if (Doc_timings)
1046 {
1047 t_send_info = TimingHelpers::timer();
1048 t_send_info_max =
1049 std::max(t_send_info_max, t_send_info - t_missing);
1050 t_send_info_min =
1051 std::min(t_send_info_min, t_send_info - t_missing);
1052 t_send_info_tot += (t_send_info - t_missing);
1053 }
1054
1055 // Create any located external halo elements on the current process
1056 create_external_halo_elements<EXT_ELEMENT>(
1057 iproc, mesh_pt, external_mesh_pt, problem_pt, interaction_index);
1058
1059 if (Doc_timings)
1060 {
1061 t_create_halo = TimingHelpers::timer();
1062 t_create_halo_max =
1063 std::max(t_create_halo_max, t_create_halo - t_send_info);
1064 t_create_halo_min =
1065 std::min(t_create_halo_min, t_create_halo - t_send_info);
1066 t_create_halo_tot += (t_create_halo - t_send_info);
1067 }
1068
1069 // Do we have any further locating to do or have we found
1070 // everything at this level of the ring communication?
1071 // Only perform the reduction operation if there's more than
1072 // one process [Note: Array is padded
1073 // by DIM times DBL_MAX entries for each mesh]
1074 n_zeta_not_found =
1075 Flat_packed_zetas_not_found_locally.size() - Dim * n_mesh;
1076
1077
1078#ifdef OOMPH_HAS_MPI
1079 if (problem_pt->communicator_pt()->nproc() > 1)
1080 {
1081 unsigned count_local_zetas = n_zeta_not_found;
1082 MPI_Allreduce(&count_local_zetas,
1083 &n_zeta_not_found,
1084 1,
1085 MPI_UNSIGNED,
1086 MPI_SUM,
1087 problem_pt->communicator_pt()->mpi_comm());
1088 }
1089#endif
1090
1091 // If its is now zero then break out of the ring comms loop
1092 if (n_zeta_not_found == 0)
1093 {
1094 if (Doc_timings)
1095 {
1096 t_local = TimingHelpers::timer();
1097 oomph_info << "BREAK N-1: CPU for entrire spiral [spiral level "
1098 << i_level << "]: " << t_local - t_spiral_start
1099 << std::endl;
1100 }
1101 break;
1102 }
1103 }
1104
1105
1106 // Doc timings
1107 if (Doc_timings)
1108 {
1109 oomph_info << "Ring-based search continued until iteration "
1110 << ring_count << " out of a maximum of "
1111 << problem_pt->communicator_pt()->nproc() - 1 << "\n";
1112 oomph_info << "Total, av, max, min CPU for send/recv of remaining "
1113 "zeta coordinates: "
1114 << t_sendrecv_tot << " "
1115 << t_sendrecv_tot / double(ring_count) << " "
1116 << t_sendrecv_max << " " << t_sendrecv_min << "\n";
1117 oomph_info << "Total, av, max, min CPU for location of missing zeta "
1118 "coordinates : "
1119 << t_missing_tot << " "
1120 << t_missing_tot / double(ring_count) << " "
1121 << t_missing_max << " " << t_missing_min << "\n";
1122 oomph_info << "Total, av, max, min CPU for send/recv of new element "
1123 "info : "
1124 << t_send_info_tot << " "
1125 << t_send_info_tot / double(ring_count) << " "
1126 << t_send_info_max << " " << t_send_info_min << "\n";
1127 oomph_info << "Total, av, max, min CPU for local creation of "
1128 "external halo objects: "
1129 << t_create_halo_tot << " "
1130 << t_create_halo_tot / double(ring_count) << " "
1131 << t_create_halo_max << " " << t_create_halo_min << "\n";
1132 }
1133
1134 } // end if for missing zetas on any process
1135#endif
1136
1137
1138 // Store information about location of elements for integration points
1139 if (Doc_stats)
1140 {
1141 unsigned count_locates = 0;
1142 unsigned n_ext_loc = External_element_located.size();
1143 for (unsigned e = 0; e < n_ext_loc; e++)
1144 {
1145 unsigned n_intpt = External_element_located[e].size();
1146 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1147 {
1148 count_locates += External_element_located[e][ipt];
1149 }
1150 }
1151
1152
1153 // Store total percentage of locates so far.
1154 // Only assign if we had anything to allocte, otherwise 100%
1155 // (default assignment) is correct
1156 if (tot_int != 0)
1157 {
1158 percentage_coords_located_elsewhere.push_back(
1159 100.0 * double(count_locates) / double(tot_int));
1160 }
1161 else
1162 {
1163 // Had none to find so we found them all!
1164 percentage_coords_located_locally.push_back(100.0);
1165 }
1166 }
1167
1168 // Do we have any further locating to do? If so, the remaining
1169 // zetas will (hopefully) be found at the next spiral level.
1170 // Only perform the reduction operation if there's more than one process
1171 // [Note: Array is padded
1172 // by DIM times DBL_MAX entries for each mesh]
1173 n_zeta_not_found =
1174 Flat_packed_zetas_not_found_locally.size() - Dim * n_mesh;
1175
1176
1177#ifdef OOMPH_HAS_MPI
1178 if (problem_pt->communicator_pt()->nproc() > 1)
1179 {
1180 unsigned count_local_zetas = n_zeta_not_found;
1181 MPI_Allreduce(&count_local_zetas,
1182 &n_zeta_not_found,
1183 1,
1184 MPI_UNSIGNED,
1185 MPI_SUM,
1186 problem_pt->communicator_pt()->mpi_comm());
1187 }
1188
1189 // Specify max level reached for later loop
1190 max_level_reached = i_level + 1;
1191#endif
1192
1193 /// If it's is now zero then break out of the spirals loop
1194 if (n_zeta_not_found == 0)
1195 {
1196 if (Doc_timings)
1197 {
1198 t_local = TimingHelpers::timer();
1199 oomph_info << "BREAK N: CPU for entrire spiral [spiral level "
1200 << i_level << "]: " << t_local - t_spiral_start
1201 << std::endl;
1202 }
1203 break;
1204 }
1205
1206 if (Doc_timings)
1207 {
1208 t_local = TimingHelpers::timer();
1209 oomph_info << "CPU for entrire spiral [spiral level " << i_level
1210 << "]: " << t_local - t_spiral_start << std::endl;
1211 }
1212
1213 // Bump up spiral levels for all meshes
1214 i_level++;
1215 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
1216 {
1217 if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
1219 {
1220 RefineableBinArray* bin_array_pt = dynamic_cast<RefineableBinArray*>(
1221 mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
1222 bin_array_pt
1224 bin_array_pt
1226 bin_array_pt
1228 bin_array_pt
1230 }
1231 else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
1233 {
1234 NonRefineableBinArray* bin_array_pt =
1235 dynamic_cast<NonRefineableBinArray*>(
1236 mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
1237
1238 bin_array_pt->current_min_spiral_level() +=
1239 bin_array_pt->n_spiral_chunk();
1240 bin_array_pt->current_max_spiral_level() +=
1241 bin_array_pt->n_spiral_chunk();
1242 }
1243#ifdef OOMPH_HAS_CGAL
1244 else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
1246 {
1247 CGALSamplePointContainer* bin_array_pt =
1248 dynamic_cast<CGALSamplePointContainer*>(
1249 mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
1250 bin_array_pt
1252 bin_array_pt
1254 bin_array_pt
1256 bin_array_pt
1258 }
1259#endif // cgal
1260 }
1261
1262 // Check termination criterion for while loop
1263 if (mesh_geom_obj_pt[0]->sample_point_container_version() ==
1265 {
1266 RefineableBinArray* bin_array_pt = dynamic_cast<RefineableBinArray*>(
1267 mesh_geom_obj_pt[0]->sample_point_container_pt());
1268
1269 if (bin_array_pt
1270 ->first_sample_point_to_actually_lookup_during_locate_zeta() <=
1271 max_n_sample_points_of_sample_point_containers)
1272 {
1273 has_not_reached_max_level_of_search = true;
1274 }
1275 else
1276 {
1277 has_not_reached_max_level_of_search = false;
1278 }
1279 }
1280 else if (mesh_geom_obj_pt[0]->sample_point_container_version() ==
1282 {
1283 NonRefineableBinArray* bin_array_pt =
1284 dynamic_cast<NonRefineableBinArray*>(
1285 mesh_geom_obj_pt[0]->sample_point_container_pt());
1286
1287 if (bin_array_pt->current_max_spiral_level() < n_max_level)
1288 {
1289 has_not_reached_max_level_of_search = true;
1290 }
1291 else
1292 {
1293 has_not_reached_max_level_of_search = false;
1294 }
1295 }
1296#ifdef OOMPH_HAS_CGAL
1297 else if (mesh_geom_obj_pt[0]->sample_point_container_version() ==
1299 {
1300 CGALSamplePointContainer* bin_array_pt =
1301 dynamic_cast<CGALSamplePointContainer*>(
1302 mesh_geom_obj_pt[0]->sample_point_container_pt());
1303
1304 if (bin_array_pt
1305 ->first_sample_point_to_actually_lookup_during_locate_zeta() <=
1306 max_n_sample_points_of_sample_point_containers)
1307 {
1308 has_not_reached_max_level_of_search = true;
1309 }
1310 else
1311 {
1312 has_not_reached_max_level_of_search = false;
1313 }
1314 }
1315#endif // cgal
1316 } // end of "spirals" loop
1317
1318
1319 // If we haven't found all zetas we're dead now...
1320 //-------------------------------------------------
1321 if (n_zeta_not_found != 0)
1322 {
1323 // Shout?
1325 {
1326 std::ostringstream error_stream;
1327 error_stream
1328 << "Multi_domain_functions::locate_zeta_for_local_coordinates()"
1329 << "\nhas failed ";
1330
1331#ifdef OOMPH_HAS_MPI
1332 if (problem_pt->communicator_pt()->nproc() > 1)
1333 {
1334 error_stream << " on proc: "
1335 << problem_pt->communicator_pt()->my_rank() << std::endl;
1336 }
1337#endif
1338 error_stream
1339 << "\n\n\nThis is most likely to arise because the two meshes\n"
1340 << "that are to be matched don't overlap perfectly or\n"
1341 << "because the elements are distorted and too small a \n"
1342 << "number of sampling points has been used when setting\n"
1343 << "up the bin structure.\n\n"
1344 << "You should try to increase the value of \n"
1345 << "the number of sample points defined in \n\n"
1346 << " "
1347 "SamplePointContainerParameters::Default_nsample_points_generated_"
1348 "per_element"
1349 << "\n\n from its current value of "
1350 << SamplePointContainerParameters::
1351 Default_nsample_points_generated_per_element
1352 << "\n"
1353 << "\n\n"
1354 << "NOTE: You can suppress this error and \"accept failure\""
1355 << " by setting the public boolean \n\n"
1356 << " "
1357 "Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_"
1358 "domain_interaction\n\n"
1359 << " to true. In this case, the pointers to external elements\n"
1360 << " that couldn't be located will remain null\n";
1361
1362 std::ostringstream modifier;
1363#ifdef OOMPH_HAS_MPI
1364 if (problem_pt->communicator_pt()->nproc() > 1)
1365 {
1366 modifier << "_proc" << problem_pt->communicator_pt()->my_rank();
1367 }
1368#endif
1369
1370 // Loop over all meshes
1371 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
1372 {
1373 // Add yet another modifier to distinguish meshes if reqd
1374 if (n_mesh > 1)
1375 {
1376 modifier << "_mesh" << i_mesh;
1377 }
1378
1379 std::ofstream outfile;
1380 char filename[100];
1381 sprintf(
1382 filename, "missing_coords_mesh%s.dat", modifier.str().c_str());
1383 outfile.open(filename);
1384 unsigned nel = mesh_pt[i_mesh]->nelement();
1385 for (unsigned e = 0; e < nel; e++)
1386 {
1387 mesh_pt[i_mesh]->finite_element_pt(e)->output(outfile);
1388 // FiniteElement::output(outfile);
1389 }
1390 outfile.close();
1391
1392 sprintf(
1393 filename, "missing_coords_ext_mesh%s.dat", modifier.str().c_str());
1394 outfile.open(filename);
1395 nel = external_mesh_pt->nelement();
1396 for (unsigned e = 0; e < nel; e++)
1397 {
1398 external_mesh_pt->finite_element_pt(e)->output(outfile);
1399 // FiniteElement::output(outfile);
1400 }
1401 outfile.close();
1402
1403 BinArray* bin_array_pt = dynamic_cast<BinArray*>(
1404 mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
1405 if (bin_array_pt != 0)
1406 {
1407 sprintf(
1408 filename, "missing_coords_bin%s.dat", modifier.str().c_str());
1409 outfile.open(filename);
1410 bin_array_pt->output_bins(outfile);
1411 outfile.close();
1412 }
1413
1414 sprintf(filename, "missing_coords%s.dat", modifier.str().c_str());
1415 outfile.open(filename);
1416 unsigned n = External_element_located.size();
1417 error_stream << "Number of unlocated elements " << n << std::endl;
1418 for (unsigned e = 0; e < n; e++)
1419 {
1420 unsigned n_intpt = External_element_located[e].size();
1421 if (n_intpt == 0)
1422 {
1423 error_stream << "Failure to locate in halo element! "
1424 << "Why are we searching there?" << std::endl;
1425 }
1426 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1427 {
1428 if (External_element_located[e][ipt] == 0)
1429 {
1430 error_stream << "Failure at element/intpt: " << e << " " << ipt
1431 << std::endl;
1432
1433 // Cast
1435 dynamic_cast<ElementWithExternalElement*>(
1436 mesh_pt[i_mesh]->element_pt(e));
1437
1438 unsigned n_dim_el = el_pt->dim();
1439 Vector<double> s(n_dim_el);
1440 for (unsigned i = 0; i < n_dim_el; i++)
1441 {
1442 s[i] = el_pt->integral_pt()->knot(ipt, i);
1443 }
1444 unsigned n_dim = el_pt->node_pt(0)->ndim();
1445 Vector<double> x(n_dim);
1446 el_pt->interpolated_x(s, x);
1447 for (unsigned i = 0; i < n_dim; i++)
1448 {
1449 error_stream << x[i] << " ";
1450 outfile << x[i] << " ";
1451 }
1452 error_stream << std::endl;
1453 outfile << std::endl;
1454 }
1455 }
1456 }
1457 outfile.close();
1458 }
1459
1460 error_stream
1461 << "Mesh and external mesh documented in missing_coords_mesh*.dat\n"
1462 << "and missing_coords_ext_mesh*.dat, respectively. Missing \n"
1463 << "coordinates in missing_coords*.dat\n";
1464 throw OomphLibError(
1465 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1466 }
1467 // Failure is deeemed to be acceptable
1468 else
1469 {
1471 << "NOTE: Haven't found " << n_zeta_not_found
1472 << " external elements in \n"
1473 << "Multi_domain_functions::aux_setup_multi_domain_interaction(...)\n"
1474 << "but this deemed to be acceptable because \n"
1475 << "Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_"
1476 "domain_interaction\n"
1477 << "is true.\n";
1478 }
1479 }
1480
1481
1482 // Doc timings if required
1483 if (Doc_timings)
1484 {
1485 t_locate = TimingHelpers::timer();
1487 << "Total CPU for location and creation of all external elements: "
1488 << t_locate - t_start << std::endl;
1489 }
1490
1491 // Delete the geometric object representing the mesh
1492 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
1493 {
1494 delete mesh_geom_obj_pt[i_mesh];
1495 }
1496
1497 // Clean up all the (extern) Vectors associated with creating the
1498 // external storage information
1499 clean_up();
1500
1501#ifdef OOMPH_HAS_MPI
1502
1503 // Output information about external storage if required
1504 if (Doc_stats)
1505 {
1506 // Report stats regarding location method
1507 bool comm_was_required = false;
1508 oomph_info << "-------------------------------------------" << std::endl;
1509 oomph_info << "- Cumulative percentage of locate success -" << std::endl;
1510 oomph_info << "--- Spiral -- Found local -- Found else ---" << std::endl;
1511 for (unsigned level = 0; level < max_level_reached; level++)
1512 {
1513 oomph_info << "--- " << level << " -- "
1514 << percentage_coords_located_locally[level] << " -- "
1515 << percentage_coords_located_elsewhere[level] << " ---"
1516 << std::endl;
1517 // Has communication with other processors at this level actually
1518 // produced any results?
1519 if (percentage_coords_located_elsewhere[level] >
1520 percentage_coords_located_locally[level])
1521 {
1522 comm_was_required = true;
1523 }
1524 }
1525 oomph_info << "-------------------------------------------" << std::endl;
1526
1527
1528 // No need for any of this malaki if we're not running in parallel
1529 if (problem_pt->communicator_pt()->nproc() > 1)
1530 {
1531 // Initialise to indicate that none of the zetas required
1532 // on this processor were located through parallel ring search,
1533 // i.e. comm was not required and we could have done some
1534 // more local searching first
1535 oomph_info << std::endl;
1536 oomph_info << "ASSESSMENT OF NEED FOR PARALLEL SEARCH: \n";
1537 oomph_info << "=======================================\n";
1538 unsigned status = 0;
1539 if (comm_was_required)
1540 {
1541 oomph_info << "- Ring-based parallel search did successfully locate "
1542 "zetas on proc "
1543 << my_rank << std::endl;
1544 status = 1;
1545 }
1546 else
1547 {
1548 if (max_level_reached > 1)
1549 {
1551 << "- Ring-based parallel search did NOT locate zetas on proc"
1552 << my_rank << std::endl;
1553 }
1554 else
1555 {
1557 << "- No ring-based parallel search was performed on proc"
1558 << my_rank << std::endl;
1559 }
1560 }
1561
1562 // Allreduce to check if anyone has benefitted from parallel ring
1563 // search
1564 unsigned overall_status = 0;
1565 MPI_Allreduce(&status,
1566 &overall_status,
1567 1,
1568 MPI_UNSIGNED,
1569 MPI_MAX,
1570 problem_pt->communicator_pt()->mpi_comm());
1571
1572 // Report of mpi was useful to anyone
1573 if (overall_status == 1)
1574 {
1575 oomph_info << "- Ring-based, parallel search did succesfully\n";
1576 oomph_info << " locate zetas on at least one other proc, so it\n";
1577 oomph_info << " was worthwhile.\n";
1578 oomph_info << std::endl;
1579 }
1580 else
1581 {
1582 if (max_level_reached > 1)
1583 {
1585 << "- Ring-based, parallel search did NOT locate zetas\n";
1586 oomph_info << " on ANY other procs, i.e it was useless.\n";
1588 << " --> We should really have done more local search\n";
1590 << " by reducing number of bins, or doing more spirals\n";
1591 oomph_info << " in one go before initiating parallel search.\n";
1592 oomph_info << std::endl;
1593 }
1594 else
1595 {
1596 oomph_info << "- No ring-based, parallel search was performed\n";
1597 oomph_info << " or necessary. Perfect!\n";
1598 oomph_info << std::endl;
1599 }
1600 }
1601 }
1602
1603 // How many external elements does the external mesh have now?
1604 oomph_info << "------------------------------------------" << std::endl;
1605 oomph_info << "External mesh: I have " << external_mesh_pt->nelement()
1606 << " elements, and " << std::endl
1607 << external_mesh_pt->nexternal_halo_element()
1608 << " external halo elements, "
1609 << external_mesh_pt->nexternal_haloed_element()
1610 << " external haloed elements" << std::endl;
1611
1612 // How many external nodes does each submesh have now?
1613 oomph_info << "------------------------------------------" << std::endl;
1614 oomph_info << "External mesh: I have " << external_mesh_pt->nnode()
1615 << " nodes, and " << std::endl
1616 << external_mesh_pt->nexternal_halo_node()
1617 << " external halo nodes, "
1618 << external_mesh_pt->nexternal_haloed_node()
1619 << " external haloed nodes" << std::endl;
1620 oomph_info << "------------------------------------------" << std::endl;
1621 }
1622
1623 // Output further information about (external) halo(ed)
1624 // elements and nodes if required
1625 if (Doc_full_stats)
1626 {
1627 // How many elements does this submesh have for each of the processors?
1628 for (int iproc = 0; iproc < n_proc; iproc++)
1629 {
1630 oomph_info << "----------------------------------------" << std::endl;
1631 oomph_info << "With process " << iproc << " there are "
1632 << external_mesh_pt->nroot_halo_element(iproc)
1633 << " root halo elements, and "
1634 << external_mesh_pt->nroot_haloed_element(iproc)
1635 << " root haloed elements" << std::endl
1636 << "and there are "
1637 << external_mesh_pt->nexternal_halo_element(iproc)
1638 << " external halo elements, and "
1639 << external_mesh_pt->nexternal_haloed_element(iproc)
1640 << " external haloed elements." << std::endl;
1641
1642 oomph_info << "----------------------------------------" << std::endl;
1643 oomph_info << "With process " << iproc << " there are "
1644 << external_mesh_pt->nhalo_node(iproc) << " halo nodes, and "
1645 << external_mesh_pt->nhaloed_node(iproc) << " haloed nodes"
1646 << std::endl
1647 << "and there are "
1648 << external_mesh_pt->nexternal_halo_node(iproc)
1649 << " external halo nodes, and "
1650 << external_mesh_pt->nexternal_haloed_node(iproc)
1651 << " external haloed nodes." << std::endl;
1652 }
1653 oomph_info << "-----------------------------------------" << std::endl
1654 << std::endl;
1655 }
1656
1657#endif
1658
1659 // Doc timings if required
1660 if (Doc_timings)
1661 {
1662 t_end = TimingHelpers::timer();
1663 oomph_info << "CPU for (one way) aux_setup_multi_domain_interaction: "
1664 << t_end - t_start << std::endl;
1665 }
1666
1667 } // end of aux_setup_multi_domain_interaction
1668
1669#ifdef OOMPH_HAS_MPI
1670
1671 //=====================================================================
1672 /// Creates external (halo) elements on the loop process based on the
1673 /// information received from each locate_zeta call on other processes.
1674 /// vector based version
1675 //=====================================================================
1676 template<class EXT_ELEMENT>
1678 int& iproc,
1679 const Vector<Mesh*>& mesh_pt,
1680 Mesh* const& external_mesh_pt,
1681 Problem* problem_pt,
1682 const unsigned& interaction_index)
1683 {
1684 OomphCommunicator* comm_pt = problem_pt->communicator_pt();
1685 int my_rank = comm_pt->my_rank();
1686
1687 // Reset counters for flat packed unsigneds (namespace data because
1688 // it's also accessed by helper functions)
1691
1692 // Initialise counter for stepping through zetas
1693 unsigned zeta_counter = 0;
1694
1695 // Initialise counter for stepping through flat-packed located
1696 // coordinates
1697 unsigned counter_for_located_coord = 0;
1698
1699 // Counter for elements in flag packed storage
1700 unsigned e_count = 0;
1701
1702 // Loop over all meshes
1703 unsigned n_mesh = mesh_pt.size();
1704 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
1705 {
1706 // The creation all happens on the current processor
1707 // Loop over this processors elements
1708 unsigned n_element = mesh_pt[i_mesh]->nelement();
1709 for (unsigned e = 0; e < n_element; e++)
1710 {
1711 // Cast to ElementWithExternalElement to set external element
1712 // (if located)
1714 dynamic_cast<ElementWithExternalElement*>(
1715 mesh_pt[i_mesh]->element_pt(e));
1716
1717 // We're not setting up external elements for halo elements
1718 // (Note: this is consistent with padding introduced when
1719 // External_element_located was first assigned)
1720 if (!el_pt->is_halo())
1721 {
1722 // Loop over integration points
1723 unsigned n_intpt = el_pt->integral_pt()->nweight();
1724 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1725 {
1726 // Has an external element been assigned to this integration point?
1727 if (External_element_located[e_count][ipt] == 0)
1728 {
1729 // Was a (non-halo) element located for this integration point
1730 if (((Proc_id_plus_one_of_external_element[zeta_counter] - 1) ==
1731 my_rank) ||
1732 (Proc_id_plus_one_of_external_element[zeta_counter] == 0))
1733 {
1734 // Either it was already found, or not found on the current
1735 // proc. In either case, we don't need to do anything for this
1736 // integration point
1737 }
1738 else
1739 {
1740 // Get the process number on which the element was located
1741 unsigned loc_p =
1742 Proc_id_plus_one_of_external_element[zeta_counter] - 1;
1743
1744 // Is it a new external halo element or not?
1745 // If so, then create it, populate it, and add it as a
1746 // source; if not, then find the right one which
1747 // has already been created and use it as the source
1748 // element.
1749
1750 // FiniteElement stored at this integration point
1751 FiniteElement* f_el_pt = 0;
1752
1753 // Is it a new element?
1754 if (Located_element_status[zeta_counter] == New)
1755 {
1756 // Create a new element from the communicated values
1757 // and coords from the process that located zeta
1758 GeneralisedElement* new_el_pt = new EXT_ELEMENT;
1759
1760 // Add external halo element to this mesh
1761 external_mesh_pt->add_external_halo_element_pt(loc_p,
1762 new_el_pt);
1763
1764 // Cast to the FE pointer
1765 f_el_pt = dynamic_cast<FiniteElement*>(new_el_pt);
1766
1767 // We need the number of interpolated values if Refineable
1768 int n_cont_inter_values = -1;
1769 if (dynamic_cast<RefineableElement*>(new_el_pt) != 0)
1770 {
1771 n_cont_inter_values =
1772 dynamic_cast<RefineableElement*>(new_el_pt)
1773 ->ncont_interpolated_values();
1774 }
1775
1776 // If we're using macro elements to update
1777#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1780 << " Using macro element node update "
1782 << std::endl;
1783#endif
1786 {
1787 // Set the macro element
1788 MacroElementNodeUpdateMesh* macro_mesh_pt =
1789 dynamic_cast<MacroElementNodeUpdateMesh*>(
1790 external_mesh_pt);
1791
1792#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1794 << " Number of macro element "
1797 << std::endl;
1798#endif
1799 unsigned macro_el_num = Flat_packed_unsigneds
1801 f_el_pt->set_macro_elem_pt(
1802 macro_mesh_pt->macro_domain_pt()->macro_element_pt(
1803 macro_el_num));
1804
1805
1806 // We need to receive the lower left
1807 // and upper right coordinates of the macro element
1808 QElementBase* q_el_pt =
1809 dynamic_cast<QElementBase*>(new_el_pt);
1810 if (q_el_pt != 0)
1811 {
1812 unsigned el_dim = q_el_pt->dim();
1813 for (unsigned i_dim = 0; i_dim < el_dim; i_dim++)
1814 {
1815 q_el_pt->s_macro_ll(i_dim) = Flat_packed_doubles
1817 q_el_pt->s_macro_ur(i_dim) = Flat_packed_doubles
1819 }
1820 }
1821 else // Throw an error, since this is only implemented for Q
1822 {
1823 std::ostringstream error_stream;
1824 error_stream << "Using MacroElement node update\n"
1825 << "in a case with non-QElements\n"
1826 << "has not yet been implemented.\n";
1827 throw OomphLibError(error_stream.str(),
1828 OOMPH_CURRENT_FUNCTION,
1829 OOMPH_EXCEPTION_LOCATION);
1830 }
1831 }
1832
1833 // Now we add nodes to the new element
1834 unsigned n_node = f_el_pt->nnode();
1835 for (unsigned j = 0; j < n_node; j++)
1836 {
1837 Node* new_nod_pt = 0;
1838
1839 // Call the add external halo node helper function
1840 add_external_halo_node_to_storage<EXT_ELEMENT>(
1841 new_nod_pt,
1842 external_mesh_pt,
1843 loc_p,
1844 j,
1845 f_el_pt,
1846 n_cont_inter_values,
1847 problem_pt);
1848 }
1849 }
1850 else // the element already exists as an external_halo
1851 {
1852#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1855 << " Index of existing external halo element "
1857 << std::endl;
1858#endif
1859 // The index itself is in Flat_packed_unsigneds[...]
1860 unsigned external_halo_el_index =
1862
1863 // Use this index to get the element
1864 f_el_pt = dynamic_cast<FiniteElement*>(
1865 external_mesh_pt->external_halo_element_pt(
1866 loc_p, external_halo_el_index));
1867
1868 // If it's not a finite element die
1869 if (f_el_pt == 0)
1870 {
1871 throw OomphLibError(
1872 "External halo element is not a FiniteElement\n",
1873 OOMPH_CURRENT_FUNCTION,
1874 OOMPH_EXCEPTION_LOCATION);
1875 }
1876 }
1877
1878 // The source element storage was initialised but
1879 // not filled earlier, so do it now
1880 // The located coordinates are required
1881 // (which could be a different dimension to zeta, e.g. in FSI)
1882 unsigned el_dim = f_el_pt->dim();
1883 Vector<double> s_located(el_dim);
1884 for (unsigned i = 0; i < el_dim; i++)
1885 {
1886 s_located[i] =
1887 Flat_packed_located_coordinates[counter_for_located_coord];
1888 counter_for_located_coord++;
1889 }
1890
1891 // Set the element for this integration point
1892 el_pt->external_element_pt(interaction_index, ipt) = f_el_pt;
1893 el_pt->external_element_local_coord(interaction_index, ipt) =
1894 s_located;
1895
1896 // Set the lookup array to true
1897 External_element_located[e_count][ipt] = 1;
1898 }
1899
1900 // Increment the integration point counter
1901 zeta_counter++;
1902 }
1903 } // end loop over integration points
1904 } // end of is halo
1905
1906 // Bump flat-packed element counter
1907 e_count++;
1908
1909 } // end of loop over elements
1910
1911 // Bump up zeta counter to skip over padding entry at end of
1912 // mesh
1913 zeta_counter++;
1914
1915 } // end loop over meshes
1916 }
1917
1918
1919 //============start of add_external_halo_node_to_storage===============
1920 /// Helper function to add external halo nodes, including any masters,
1921 /// based on information received from the haloed process
1922 //=========================================================================
1923 template<class EXT_ELEMENT>
1925 Node*& new_nod_pt,
1926 Mesh* const& external_mesh_pt,
1927 unsigned& loc_p,
1928 unsigned& node_index,
1929 FiniteElement* const& new_el_pt,
1930 int& n_cont_inter_values,
1931 Problem* problem_pt)
1932 {
1933 // Add the external halo node if required
1935 external_mesh_pt,
1936 loc_p,
1937 node_index,
1938 new_el_pt,
1939 n_cont_inter_values,
1940 problem_pt);
1941
1942 // Recursively add masters
1943 recursively_add_masters_of_external_halo_node_to_storage<EXT_ELEMENT>(
1944 new_nod_pt,
1945 external_mesh_pt,
1946 loc_p,
1947 node_index,
1948 new_el_pt,
1949 n_cont_inter_values,
1950 problem_pt);
1951 }
1952
1953
1954 //========================================================================
1955 /// Recursively add masters of external halo nodes (and their masters, etc)
1956 /// based on information received from the haloed process
1957 //=========================================================================
1958 template<class EXT_ELEMENT>
1961 Node*& new_nod_pt,
1962 Mesh* const& external_mesh_pt,
1963 unsigned& loc_p,
1964 unsigned& node_index,
1965 FiniteElement* const& new_el_pt,
1966 int& n_cont_inter_values,
1967 Problem* problem_pt)
1968 {
1969 for (int i_cont = -1; i_cont < n_cont_inter_values; i_cont++)
1970 {
1971#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1973 << " Boolean to indicate that continuously interpolated "
1974 "variable i_cont "
1975 << i_cont << " is hanging "
1977 << std::endl;
1978#endif
1980 {
1981#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1983 << " Number of master nodes "
1985 << std::endl;
1986#endif
1987 unsigned n_master =
1989
1990 // Setup new HangInfo
1991 HangInfo* hang_pt = new HangInfo(n_master);
1992 for (unsigned m = 0; m < n_master; m++)
1993 {
1994 Node* master_nod_pt = 0;
1995 // Get the master node (creating and adding it if required)
1996 add_external_halo_master_node_helper<EXT_ELEMENT>(master_nod_pt,
1997 new_nod_pt,
1998 external_mesh_pt,
1999 loc_p,
2000 n_cont_inter_values,
2001 problem_pt);
2002
2003 // Get the weight and set the HangInfo
2004 double master_weight =
2006 hang_pt->set_master_node_pt(m, master_nod_pt, master_weight);
2007
2008 // Recursively add masters of master
2009 recursively_add_masters_of_external_halo_node_to_storage<EXT_ELEMENT>(
2010 master_nod_pt,
2011 external_mesh_pt,
2012 loc_p,
2013 node_index,
2014 new_el_pt,
2015 n_cont_inter_values,
2016 problem_pt);
2017 }
2018 new_nod_pt->set_hanging_pt(hang_pt, i_cont);
2019 }
2020 } // end loop over continous interpolated values
2021 }
2022
2023 //========================================================================
2024 /// Helper function to add external halo node that is a master
2025 //========================================================================
2026 template<class EXT_ELEMENT>
2028 Node*& new_master_nod_pt,
2029 Node*& new_nod_pt,
2030 Mesh* const& external_mesh_pt,
2031 unsigned& loc_p,
2032 int& ncont_inter_values,
2033 Problem* problem_pt)
2034 {
2035 // Given the node and the external mesh, and received information
2036 // about them from process loc_p, construct them on the current process
2037#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2040 << " Boolean to trigger construction of new external halo master node "
2042#endif
2044 {
2045 // Construct a new node based upon sent information
2046 construct_new_external_halo_master_node_helper<EXT_ELEMENT>(
2047 new_master_nod_pt, new_nod_pt, loc_p, external_mesh_pt, problem_pt);
2048 }
2049 else
2050 {
2051#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2053 << " index of existing external halo master node "
2055 << std::endl;
2056#endif
2057 // Copy node from received location
2058 new_master_nod_pt = external_mesh_pt->external_halo_node_pt(
2060 }
2061 }
2062
2063 //======start of construct_new_external_halo_master_node_helper===========
2064 /// Helper function which constructs a new external halo master node
2065 /// with the required information sent from the haloed process
2066 //========================================================================
2067 template<class EXT_ELEMENT>
2069 Node*& new_master_nod_pt,
2070 Node*& nod_pt,
2071 unsigned& loc_p,
2072 Mesh* const& external_mesh_pt,
2073 Problem* problem_pt)
2074 {
2075 // First three sent numbers are dimension, position type and nvalue
2076 // (to be used in Node constructors)
2077#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2079 << " ndim for external halo master node "
2081 << std::endl;
2082#endif
2084#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2086 << " nposition type for external halo master node "
2088 << std::endl;
2089#endif
2090 unsigned n_position_type =
2092#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2094 << " nvalue for external halo master node "
2096 << std::endl;
2097#endif
2098 unsigned n_value =
2100
2101 // If it's a solid node also receive the lagrangian dimension and pos type
2102 SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
2103 unsigned n_lag_dim;
2104 unsigned n_lag_type;
2105 if (solid_nod_pt != 0)
2106 {
2107#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2109 << " nlagrdim for external halo master solid node "
2111 << std::endl;
2112#endif
2114#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2116 << " nlagrtype for external halo master solid node "
2118 << std::endl;
2119#endif
2121 }
2122
2123 // Null TimeStepper for now
2124 TimeStepper* time_stepper_pt = 0;
2125 // Default number of previous values to 1
2126 unsigned n_prev = 1;
2127
2128 // The first entry in nodal_info indicates
2129 // the timestepper required for this halo node
2130#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2132 << " Boolean: need timestepper "
2134 << std::endl;
2135#endif
2137 {
2138#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2140 << " Index minus one of timestepper "
2142 << std::endl;
2143#endif
2144 // Index minus one!
2145 time_stepper_pt = problem_pt->time_stepper_pt(
2147
2148 // Check whether number of prev values is "sent" across
2149 n_prev = time_stepper_pt->ntstorage();
2150 }
2151
2152 // Is the node for which the master is required Algebraic, Macro or Solid?
2153 AlgebraicNode* alg_nod_pt = dynamic_cast<AlgebraicNode*>(nod_pt);
2154 MacroElementNodeUpdateNode* macro_nod_pt =
2155 dynamic_cast<MacroElementNodeUpdateNode*>(nod_pt);
2156
2157 // What type of node was the node for which we are constructing a master?
2158 if (alg_nod_pt != 0)
2159 {
2160 // The master node should also be algebraic
2161#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2163 << " Boolean for algebraic boundary node "
2165 << std::endl;
2166#endif
2167 // If this master node's haloed copy is on a boundary then
2168 // it needs to be on the same boundary here
2170 {
2171 // Create a new BoundaryNode (not attached to an element)
2172 if (time_stepper_pt != 0)
2173 {
2174 new_master_nod_pt = new BoundaryNode<AlgebraicNode>(
2175 time_stepper_pt, n_dim, n_position_type, n_value);
2176 }
2177 else
2178 {
2179 new_master_nod_pt =
2180 new BoundaryNode<AlgebraicNode>(n_dim, n_position_type, n_value);
2181 }
2182
2183 // How many boundaries does the algebraic master node live on?
2184#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2186 << " Number of boundaries the algebraic master node is on: "
2188 << std::endl;
2189#endif
2190 unsigned nb =
2192 for (unsigned i = 0; i < nb; i++)
2193 {
2194 // Boundary number
2195#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2197 << " Algebraic master node is on boundary "
2199 << std::endl;
2200#endif
2201 unsigned i_bnd =
2203 external_mesh_pt->add_boundary_node(i_bnd, new_master_nod_pt);
2204 }
2205
2206
2207 // Do we have additional values created by face elements?
2208#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2210 << "Number of additional values created by face element "
2211 << "for master node "
2213 << std::endl;
2214#endif
2215 unsigned n_entry =
2217 if (n_entry > 0)
2218 {
2219 // Create storage, if it doesn't already exist, for the map
2220 // that will contain the position of the first entry of
2221 // this face element's additional values,
2222 BoundaryNodeBase* bnew_master_nod_pt =
2223 dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
2224#ifdef PARANOID
2225 if (bnew_master_nod_pt == 0)
2226 {
2227 throw OomphLibError("Failed to cast new node to boundary node\n",
2228 OOMPH_CURRENT_FUNCTION,
2229 OOMPH_EXCEPTION_LOCATION);
2230 }
2231#endif
2232 if (bnew_master_nod_pt
2233 ->index_of_first_value_assigned_by_face_element_pt() == 0)
2234 {
2235 bnew_master_nod_pt
2237 new std::map<unsigned, unsigned>;
2238 }
2239
2240 // Get pointer to the map of indices associated with
2241 // additional values created by face elements
2242 std::map<unsigned, unsigned>* map_pt =
2243 bnew_master_nod_pt
2245
2246 // Loop over number of entries in map
2247 for (unsigned i = 0; i < n_entry; i++)
2248 {
2249 // Read out pairs...
2250
2251#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2254 << " Key of map entry for master node"
2256 << std::endl;
2257#endif
2258 unsigned first =
2260
2261#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2264 << " Value of map entry for master node"
2266 << std::endl;
2267#endif
2268 unsigned second =
2270
2271 // ...and assign
2272 (*map_pt)[first] = second;
2273 }
2274 }
2275 }
2276 else
2277 {
2278 // Create node (not attached to any element)
2279 if (time_stepper_pt != 0)
2280 {
2281 new_master_nod_pt =
2282 new AlgebraicNode(time_stepper_pt, n_dim, n_position_type, n_value);
2283 }
2284 else
2285 {
2286 new_master_nod_pt =
2287 new AlgebraicNode(n_dim, n_position_type, n_value);
2288 }
2289 }
2290
2291 // Add this as an external halo node BEFORE considering node update!
2292 external_mesh_pt->add_external_halo_node_pt(loc_p, new_master_nod_pt);
2293
2294 // The external mesh is itself Algebraic...
2295 AlgebraicMesh* alg_mesh_pt =
2296 dynamic_cast<AlgebraicMesh*>(external_mesh_pt);
2297
2298#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2300 << " algebraic node update id for master node "
2302 << std::endl;
2303#endif
2304 /// The first entry of All_unsigned_values is the default node update id
2305 unsigned update_id =
2307
2308 // Setup algebraic node update info for this new node
2309 Vector<double> ref_value;
2310
2311#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2313 << " algebraic node number of ref values for master node "
2315 << std::endl;
2316#endif
2317 // The size of this vector is in the next entry
2318 unsigned n_ref_val =
2320
2321 // The reference values are in the subsequent entries of All_double_values
2322 ref_value.resize(n_ref_val);
2323 for (unsigned i_ref = 0; i_ref < n_ref_val; i_ref++)
2324 {
2325 ref_value[i_ref] =
2327 }
2328
2329 // Also require a Vector of geometric objects
2330 Vector<GeomObject*> geom_object_pt;
2331
2332#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2334 << " algebraic node number of geom objects for master node "
2336 << std::endl;
2337#endif
2338
2339 // The size of this vector is in the next entry of All_unsigned_values
2340 unsigned n_geom_obj =
2342
2343 // The remaining indices are in the rest of
2344 // All_alg_nodal_info
2345 geom_object_pt.resize(n_geom_obj);
2346 for (unsigned i_geom = 0; i_geom < n_geom_obj; i_geom++)
2347 {
2348#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2350 << " algebraic node: " << i_geom << "-th out of "
2351 << n_geom_obj << "-th geom index "
2353 << std::endl;
2354#endif
2355 unsigned geom_index =
2357
2358 // This index indicates which (if any) of the AlgebraicMesh's
2359 // stored geometric objects should be used
2360 geom_object_pt[i_geom] = alg_mesh_pt->geom_object_list_pt(geom_index);
2361 }
2362
2363 AlgebraicNode* alg_master_nod_pt =
2364 dynamic_cast<AlgebraicNode*>(new_master_nod_pt);
2365
2366 /// ... so for the specified update_id, call
2367 /// add_node_update_info
2368 alg_master_nod_pt->add_node_update_info(
2369 update_id, alg_mesh_pt, geom_object_pt, ref_value);
2370
2371 /// Now call update_node_update
2372 alg_mesh_pt->update_node_update(alg_master_nod_pt);
2373 }
2374 else if (macro_nod_pt != 0)
2375 {
2376 // The master node should also be a macro node
2377 // If this master node's haloed copy is on a boundary then
2378 // it needs to be on the same boundary here
2379#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2381 << " Boolean for master algebraic node is boundary node "
2383 << std::endl;
2384#endif
2386 {
2387 // Create a new BoundaryNode (not attached to an element)
2388 if (time_stepper_pt != 0)
2389 {
2390 new_master_nod_pt = new BoundaryNode<MacroElementNodeUpdateNode>(
2391 time_stepper_pt, n_dim, n_position_type, n_value);
2392 }
2393 else
2394 {
2395 new_master_nod_pt = new BoundaryNode<MacroElementNodeUpdateNode>(
2396 n_dim, n_position_type, n_value);
2397 }
2398
2399
2400 // How many boundaries does the macro element master node live on?
2401#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2404 << " Number of boundaries the macro element master node is on: "
2406 << std::endl;
2407#endif
2408 unsigned nb =
2410 for (unsigned i = 0; i < nb; i++)
2411 {
2412 // Boundary number
2413#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2415 << " Macro element master node is on boundary "
2417 << std::endl;
2418#endif
2419 unsigned i_bnd =
2421 external_mesh_pt->add_boundary_node(i_bnd, new_master_nod_pt);
2422 }
2423
2424 // Do we have additional values created by face elements?
2425#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2427 << " Number of additional values created by face element "
2428 << "for macro element master node "
2430 << std::endl;
2431#endif
2432 unsigned n_entry =
2434 if (n_entry > 0)
2435 {
2436 // Create storage, if it doesn't already exist, for the map
2437 // that will contain the position of the first entry of
2438 // this face element's additional values,
2439 BoundaryNodeBase* bnew_master_nod_pt =
2440 dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
2441#ifdef PARANOID
2442 if (bnew_master_nod_pt == 0)
2443 {
2444 throw OomphLibError("Failed to cast new node to boundary node\n",
2445 OOMPH_CURRENT_FUNCTION,
2446 OOMPH_EXCEPTION_LOCATION);
2447 }
2448#endif
2449 if (bnew_master_nod_pt
2450 ->index_of_first_value_assigned_by_face_element_pt() == 0)
2451 {
2452 bnew_master_nod_pt
2454 new std::map<unsigned, unsigned>;
2455 }
2456
2457 // Get pointer to the map of indices associated with
2458 // additional values created by face elements
2459 std::map<unsigned, unsigned>* map_pt =
2460 bnew_master_nod_pt
2462
2463 // Loop over number of entries in map
2464 for (unsigned i = 0; i < n_entry; i++)
2465 {
2466 // Read out pairs...
2467
2468#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2471 << " Key of map entry for macro element master node"
2473 << std::endl;
2474#endif
2475 unsigned first =
2477
2478#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2481 << " Value of map entry for macro element master node"
2483 << std::endl;
2484#endif
2485 unsigned second =
2487
2488 // ...and assign
2489 (*map_pt)[first] = second;
2490 }
2491 }
2492 }
2493 else
2494 {
2495 // Create node (not attached to any element)
2496 if (time_stepper_pt != 0)
2497 {
2498 new_master_nod_pt = new MacroElementNodeUpdateNode(
2499 time_stepper_pt, n_dim, n_position_type, n_value);
2500 }
2501 else
2502 {
2503 new_master_nod_pt =
2504 new MacroElementNodeUpdateNode(n_dim, n_position_type, n_value);
2505 }
2506 }
2507
2508 // Add this as an external halo node
2509 external_mesh_pt->add_external_halo_node_pt(loc_p, new_master_nod_pt);
2510
2511 // Create a new node update element for this master node if required
2512 FiniteElement* new_node_update_f_el_pt = 0;
2513#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2515 << " Bool: need new external halo element "
2517 << std::endl;
2518#endif
2520 {
2521 GeneralisedElement* new_node_update_el_pt = new EXT_ELEMENT;
2522
2523 // Add external hal element to this mesh
2524 external_mesh_pt->add_external_halo_element_pt(loc_p,
2525 new_node_update_el_pt);
2526
2527 // Cast to finite element
2528 new_node_update_f_el_pt =
2529 dynamic_cast<FiniteElement*>(new_node_update_el_pt);
2530
2531 // Need number of interpolated values if Refineable
2532 int n_cont_inter_values;
2533 if (dynamic_cast<RefineableElement*>(new_node_update_f_el_pt) != 0)
2534 {
2535 n_cont_inter_values =
2536 dynamic_cast<RefineableElement*>(new_node_update_f_el_pt)
2537 ->ncont_interpolated_values();
2538 }
2539 else
2540 {
2541 n_cont_inter_values = -1;
2542 }
2543#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2545 << " Bool: we have a macro element mesh "
2547 << std::endl;
2548#endif
2549 // If we're using macro elements to update,
2551 {
2552 // Set the macro element
2553 MacroElementNodeUpdateMesh* macro_mesh_pt =
2554 dynamic_cast<MacroElementNodeUpdateMesh*>(external_mesh_pt);
2555
2556#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2558 << " Number of macro element "
2560 << std::endl;
2561#endif
2562 unsigned macro_el_num =
2564 new_node_update_f_el_pt->set_macro_elem_pt(
2565 macro_mesh_pt->macro_domain_pt()->macro_element_pt(macro_el_num));
2566
2567 // we need to receive
2568 // the lower left and upper right coordinates of the macro
2569 QElementBase* q_el_pt =
2570 dynamic_cast<QElementBase*>(new_node_update_f_el_pt);
2571 if (q_el_pt != 0)
2572 {
2573 unsigned el_dim = q_el_pt->dim();
2574 for (unsigned i_dim = 0; i_dim < el_dim; i_dim++)
2575 {
2576 q_el_pt->s_macro_ll(i_dim) =
2578 q_el_pt->s_macro_ur(i_dim) =
2580 }
2581 }
2582 else // Throw an error
2583 {
2584 std::ostringstream error_stream;
2585 error_stream << "You are using a MacroElement node update\n"
2586 << "in a case with non-QElements. This has not\n"
2587 << "yet been implemented.\n";
2588 throw OomphLibError(error_stream.str(),
2589 OOMPH_CURRENT_FUNCTION,
2590 OOMPH_EXCEPTION_LOCATION);
2591 }
2592 }
2593
2594 unsigned n_node = new_node_update_f_el_pt->nnode();
2595 for (unsigned j = 0; j < n_node; j++)
2596 {
2597 Node* new_nod_pt = 0;
2598 add_external_halo_node_to_storage<EXT_ELEMENT>(
2599 new_nod_pt,
2600 external_mesh_pt,
2601 loc_p,
2602 j,
2603 new_node_update_f_el_pt,
2604 n_cont_inter_values,
2605 problem_pt);
2606 }
2607 }
2608 else // The node update element exists already
2609 {
2610#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2612 << " Number of already existing external halo element "
2614 << std::endl;
2615#endif
2616 new_node_update_f_el_pt = dynamic_cast<FiniteElement*>(
2617 external_mesh_pt->external_halo_element_pt(
2619 }
2620
2621 // Remaining required information to create functioning
2622 // MacroElementNodeUpdateNode...
2623
2624 // Get the required geom objects for the node update
2625 // from the mesh
2626 Vector<GeomObject*> geom_object_vector_pt;
2627 MacroElementNodeUpdateMesh* macro_mesh_pt =
2628 dynamic_cast<MacroElementNodeUpdateMesh*>(external_mesh_pt);
2629 geom_object_vector_pt = macro_mesh_pt->geom_object_vector_pt();
2630
2631 // Cast to MacroElementNodeUpdateNode
2632 MacroElementNodeUpdateNode* macro_master_nod_pt =
2633 dynamic_cast<MacroElementNodeUpdateNode*>(new_master_nod_pt);
2634
2635 // Set all required information - node update element,
2636 // local coordinate in this element, and then set node update info
2637 macro_master_nod_pt->node_update_element_pt() = new_node_update_f_el_pt;
2638
2639 // Need to get the local node index of the macro_master_nod_pt
2640 unsigned local_node_index;
2641 unsigned n_node = new_node_update_f_el_pt->nnode();
2642 for (unsigned j = 0; j < n_node; j++)
2643 {
2644 if (macro_master_nod_pt == new_node_update_f_el_pt->node_pt(j))
2645 {
2646 local_node_index = j;
2647 break;
2648 }
2649 }
2650
2651 Vector<double> s_in_macro_node_update_element;
2652 new_node_update_f_el_pt->local_coordinate_of_node(
2653 local_node_index, s_in_macro_node_update_element);
2654
2655 macro_master_nod_pt->set_node_update_info(new_node_update_f_el_pt,
2656 s_in_macro_node_update_element,
2657 geom_object_vector_pt);
2658 }
2659 else if (solid_nod_pt != 0)
2660 {
2661 // The master node should also be a SolidNode
2662 // If this node was on a boundary then it needs to
2663 // be on the same boundary here
2664#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2666 << " Bool master is a boundary (solid) node "
2668 << std::endl;
2669#endif
2671 {
2672 // Construct a new boundary node
2673 if (time_stepper_pt != 0)
2674 {
2675 new_master_nod_pt = new BoundaryNode<SolidNode>(time_stepper_pt,
2676 n_lag_dim,
2677 n_lag_type,
2678 n_dim,
2679 n_position_type,
2680 n_value);
2681 }
2682 else
2683 {
2684 new_master_nod_pt = new BoundaryNode<SolidNode>(
2685 n_lag_dim, n_lag_type, n_dim, n_position_type, n_value);
2686 }
2687
2688
2689 // How many boundaries does the macro element master node live on?
2690#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2692 << " Number of boundaries the solid master node is on: "
2694 << std::endl;
2695#endif
2696 unsigned nb =
2698 for (unsigned i = 0; i < nb; i++)
2699 {
2700 // Boundary number
2701#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2703 << " Solid master node is on boundary "
2705 << std::endl;
2706#endif
2707 unsigned i_bnd =
2709 external_mesh_pt->add_boundary_node(i_bnd, new_master_nod_pt);
2710 }
2711
2712 // Do we have additional values created by face elements?
2713#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2715 << " Number of additional values created by face element "
2716 << "for solid master node "
2718 << std::endl;
2719#endif
2720 unsigned n_entry =
2722 if (n_entry > 0)
2723 {
2724 // Create storage, if it doesn't already exist, for the map
2725 // that will contain the position of the first entry of
2726 // this face element's additional values,
2727 BoundaryNodeBase* bnew_master_nod_pt =
2728 dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
2729#ifdef PARANOID
2730 if (bnew_master_nod_pt == 0)
2731 {
2732 throw OomphLibError("Failed to cast new node to boundary node\n",
2733 OOMPH_CURRENT_FUNCTION,
2734 OOMPH_EXCEPTION_LOCATION);
2735 }
2736#endif
2737 if (bnew_master_nod_pt
2738 ->index_of_first_value_assigned_by_face_element_pt() == 0)
2739 {
2740 bnew_master_nod_pt
2742 new std::map<unsigned, unsigned>;
2743 }
2744
2745 // Get pointer to the map of indices associated with
2746 // additional values created by face elements
2747 std::map<unsigned, unsigned>* map_pt =
2748 bnew_master_nod_pt
2750
2751 // Loop over number of entries in map
2752 for (unsigned i = 0; i < n_entry; i++)
2753 {
2754 // Read out pairs...
2755
2756#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2759 << " Key of map entry for solid master node"
2761 << std::endl;
2762#endif
2763 unsigned first =
2765
2766#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2769 << " Value of map entry for solid master node"
2771 << std::endl;
2772#endif
2773 unsigned second =
2775
2776 // ...and assign
2777 (*map_pt)[first] = second;
2778 }
2779 }
2780 }
2781 else
2782 {
2783 // Construct an ordinary (non-boundary) node
2784 if (time_stepper_pt != 0)
2785 {
2786 new_master_nod_pt = new SolidNode(time_stepper_pt,
2787 n_lag_dim,
2788 n_lag_type,
2789 n_dim,
2790 n_position_type,
2791 n_value);
2792 }
2793 else
2794 {
2795 new_master_nod_pt = new SolidNode(
2796 n_lag_dim, n_lag_type, n_dim, n_position_type, n_value);
2797 }
2798 }
2799
2800 // Add this as an external halo node
2801 external_mesh_pt->add_external_halo_node_pt(loc_p, new_master_nod_pt);
2802
2803 // Copy across particular info required for SolidNode
2804 // NOTE: Are there any problems with additional values for SolidNodes?
2805 SolidNode* solid_master_nod_pt =
2806 dynamic_cast<SolidNode*>(new_master_nod_pt);
2807 unsigned n_solid_val =
2808 solid_master_nod_pt->variable_position_pt()->nvalue();
2809 for (unsigned i_val = 0; i_val < n_solid_val; i_val++)
2810 {
2811 for (unsigned t = 0; t < n_prev; t++)
2812 {
2813 solid_master_nod_pt->variable_position_pt()->set_value(
2815 }
2816 }
2817 }
2818 else // Just an ordinary node!
2819 {
2820#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2822 << " Bool node is on boundary "
2824 << std::endl;
2825#endif
2826
2827 // If this node was on a boundary then it needs to
2828 // be on the same boundary here
2830 {
2831 // Construct a new boundary node
2832 if (time_stepper_pt != 0)
2833 {
2834 new_master_nod_pt = new BoundaryNode<Node>(
2835 time_stepper_pt, n_dim, n_position_type, n_value);
2836 }
2837 else
2838 {
2839 new_master_nod_pt =
2840 new BoundaryNode<Node>(n_dim, n_position_type, n_value);
2841 }
2842
2843 // How many boundaries does the master node live on?
2844#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2846 << " Number of boundaries the master node is on: "
2848 << std::endl;
2849#endif
2850 unsigned nb =
2852 for (unsigned i = 0; i < nb; i++)
2853 {
2854 // Boundary number
2855#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2857 << " Master node is on boundary "
2859 << std::endl;
2860#endif
2861 unsigned i_bnd =
2863 external_mesh_pt->add_boundary_node(i_bnd, new_master_nod_pt);
2864 }
2865
2866
2867 // Do we have additional values created by face elements?
2868#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2870 << " Number of additional values created by face element "
2871 << "for master node "
2873 << std::endl;
2874#endif
2875 unsigned n_entry =
2877 if (n_entry > 0)
2878 {
2879 // Create storage, if it doesn't already exist, for the map
2880 // that will contain the position of the first entry of
2881 // this face element's additional values,
2882 BoundaryNodeBase* bnew_master_nod_pt =
2883 dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
2884#ifdef PARANOID
2885 if (bnew_master_nod_pt == 0)
2886 {
2887 throw OomphLibError("Failed to cast new node to boundary node\n",
2888 OOMPH_CURRENT_FUNCTION,
2889 OOMPH_EXCEPTION_LOCATION);
2890 }
2891#endif
2892 if (bnew_master_nod_pt
2893 ->index_of_first_value_assigned_by_face_element_pt() == 0)
2894 {
2895 bnew_master_nod_pt
2897 new std::map<unsigned, unsigned>;
2898 }
2899
2900 // Get pointer to the map of indices associated with
2901 // additional values created by face elements
2902 std::map<unsigned, unsigned>* map_pt =
2903 bnew_master_nod_pt
2905
2906 // Loop over number of entries in map
2907 for (unsigned i = 0; i < n_entry; i++)
2908 {
2909 // Read out pairs...
2910
2911#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2914 << " Key of map entry for master node"
2916 << std::endl;
2917#endif
2918 unsigned first =
2920
2921#ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2924 << " Value of map entry for master node"
2926 << std::endl;
2927#endif
2928 unsigned second =
2930
2931 // ...and assign
2932 (*map_pt)[first] = second;
2933 }
2934 }
2935 }
2936 else
2937 {
2938 // Construct an ordinary (non-boundary) node
2939 if (time_stepper_pt != 0)
2940 {
2941 new_master_nod_pt =
2942 new Node(time_stepper_pt, n_dim, n_position_type, n_value);
2943 }
2944 else
2945 {
2946 new_master_nod_pt = new Node(n_dim, n_position_type, n_value);
2947 }
2948 }
2949
2950 // Add this as an external halo node
2951 external_mesh_pt->add_external_halo_node_pt(loc_p, new_master_nod_pt);
2952 }
2953
2954 // Remaining info received for all node types
2955 // Get copied history values
2956 // unsigned n_val=new_master_nod_pt->nvalue();
2957 for (unsigned i_val = 0; i_val < n_value; i_val++)
2958 {
2959 for (unsigned t = 0; t < n_prev; t++)
2960 {
2961 new_master_nod_pt->set_value(
2963 }
2964 }
2965
2966 // Get copied history values for positions
2967 unsigned n_nod_dim = new_master_nod_pt->ndim();
2968 for (unsigned idim = 0; idim < n_nod_dim; idim++)
2969 {
2970 for (unsigned t = 0; t < n_prev; t++)
2971 {
2972 // Copy to coordinate
2973 new_master_nod_pt->x(t, idim) =
2975 }
2976 }
2977 }
2978
2979
2980#endif
2981
2982
2983} // namespace oomph
2984
2985#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
unsigned max_bin_dimension() const
Max. bin dimension (number of bins in coordinate directions)
virtual void output_bins(std::ofstream &outfile)=0
Output bins (boundaries and number of sample points in them)
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
unsigned & multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta()
Every time we've completed a "spiral", visiting a finite number of sample points in a deterministic o...
unsigned & initial_last_sample_point_to_actually_lookup_during_locate_zeta()
When searching through sample points only actually do the locate_zeta calls when when the counter exc...
unsigned & last_sample_point_to_actually_lookup_during_locate_zeta()
When searching through sample points only actually do the locate_zeta calls when when the counter is ...
unsigned & first_sample_point_to_actually_lookup_during_locate_zeta()
When searching through sample points only actually do the locate_zeta calls when when the counter exc...
unsigned total_number_of_sample_points_computed_recursively() const
Compute total number of sample points in sample point container.
///////////////////////////////////////////////////////////////////////////// ///////////////////////...
unsigned n_spiral_chunk() const
Number of spirals to be searched in one go const version.
unsigned & current_min_spiral_level()
Access function to current min. spiral level.
unsigned & current_max_spiral_level()
Access function to current max. spiral level.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
unsigned & multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta()
Every time we've completed a "spiral", visiting a finite number of sample points in a deterministic o...
unsigned total_number_of_sample_points_computed_recursively() const
Compute total number of sample points recursively.
unsigned & last_sample_point_to_actually_lookup_during_locate_zeta()
When searching through sample points recursively from the top level RefineableBinArray (in determinis...
unsigned & initial_last_sample_point_to_actually_lookup_during_locate_zeta()
When searching through sample points recursively from the top level RefineableBinArray (in determinis...
unsigned & first_sample_point_to_actually_lookup_during_locate_zeta()
When searching through sample points recursively from the top level RefineableBinArray (in determinis...
////////////////////////////////////////////////////////////////////
virtual void update_node_update(AlgebraicNode *&node_pt)=0
Update the node update info for given node, following mesh adaptation. Must be implemented for every ...
GeomObject * geom_object_list_pt(const unsigned &i)
Access function to the ith GeomObject.
////////////////////////////////////////////////////////////////////
void add_node_update_info(const int &id, AlgebraicMesh *mesh_pt, const Vector< GeomObject * > &geom_object_pt, const Vector< double > &ref_value, const bool &called_from_constructor=false)
Add algebraic update information for node: What's the ID of the mesh update function (typically used ...
A class that contains the information required by Nodes that are located on Mesh boundaries....
Definition: nodes.h:1996
std::map< unsigned, unsigned > *& index_of_first_value_assigned_by_face_element_pt()
Return pointer to the map giving the index of the first face element value.
Definition: nodes.h:2046
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh....
Definition: nodes.h:2242
A class to do comparison of the elements by lexicographic ordering, based on the boundary coordinates...
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
MacroElement * macro_element_pt(const unsigned &i)
Access to i-th macro element.
Definition: domain.h:116
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void initialise_external_element_storage()
Initialise storage for pointers to external elements and their local coordinates. This must be called...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
Class that is used to create FaceElement from bulk elements and to provide these FaceElement with a g...
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
void set_boundary_number_in_bulk_mesh(const unsigned &b)
Set function for the boundary number in bulk mesh.
Definition: elements.h:4482
A general Finite Element class.
Definition: elements.h:1313
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Set pointer to macro element – can be overloaded in derived elements to perform additional tasks.
Definition: elements.h:1872
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
Definition: elements.h:1842
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4675
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
A Generalised Element class.
Definition: elements.h:73
bool is_halo() const
Is this element a halo?
Definition: elements.h:1163
Class that contains data for hanging nodes.
Definition: nodes.h:742
void set_master_node_pt(const unsigned &i, Node *const &master_node_pt, const double &weight)
Set the pointer to the i-th master node and its weight.
Definition: nodes.cc:1474
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
MacroElementNodeUpdateMeshes contain MacroElementNodeUpdateNodes which have their own node update fun...
Vector< GeomObject * > geom_object_vector_pt()
Access function to the vector of GeomObject.
Domain *& macro_domain_pt()
Broken assignment operator.
////////////////////////////////////////////////////////////////////
void set_node_update_info(FiniteElement *node_update_element_pt, const Vector< double > &s_in_node_update_element, const Vector< GeomObject * > &geom_object_pt)
Set node update information for node: Pass the pointer to the element that performs the update operat...
FiniteElement *& node_update_element_pt()
Pointer to finite element that performs the update by referring to its macro-element representation (...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
A general mesh class.
Definition: mesh.h:67
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
Node *& external_halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external halo node in this Mesh whose non-halo external counterpart is held on...
Definition: mesh.h:2377
GeneralisedElement *& external_halo_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th external halo element in this Mesh whose non-halo counterpart is held on proce...
Definition: mesh.h:2251
unsigned nexternal_haloed_element()
Total number of external haloed elements in this Mesh.
Definition: mesh.h:2267
unsigned nexternal_halo_node()
Total number of external halo nodes in this Mesh.
Definition: mesh.h:2308
void add_external_halo_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add external halo element whose non-halo counterpart is held on processor p to this Mesh.
Definition: mesh.h:2259
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
unsigned nroot_haloed_element()
Total number of root haloed elements in this Mesh.
Definition: mesh.h:1921
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
void delete_all_external_storage()
Wipe the storage for all externally-based elements.
Definition: mesh.cc:9190
unsigned nhalo_node()
Total number of halo nodes in this Mesh.
Definition: mesh.h:1878
void add_external_halo_node_pt(const unsigned &p, Node *&nod_pt)
Add external halo node whose non-halo (external) counterpart is held on processor p to the storage sc...
Definition: mesh.h:2368
unsigned nhaloed_node()
Total number of haloed nodes in this Mesh.
Definition: mesh.h:1985
unsigned nexternal_haloed_node()
Total number of external haloed nodes in this Mesh.
Definition: mesh.h:2412
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
unsigned nexternal_halo_element()
Total number of external halo elements in this Mesh.
Definition: mesh.h:2222
unsigned nroot_halo_element()
Total number of root halo elements in this Mesh.
Definition: mesh.h:1830
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
void set_hanging_pt(HangInfo *const &hang_pt, const int &i)
Set the hanging data for the i-th value. (hang_pt=0 to make non-hanging)
Definition: nodes.cc:2068
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1246
bool problem_has_been_distributed()
Access to Problem_has_been_distributed flag.
Definition: problem.h:2178
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1524
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition: Qelements.h:91
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:202
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:186
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1765
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
bool Doc_full_stats
Boolean to indicate whether to output further info during setup_multi_domain_interaction() routines.
void add_external_halo_node_helper(Node *&new_nod_pt, Mesh *const &mesh_pt, unsigned &loc_p, unsigned &node_index, FiniteElement *const &new_el_pt, int &n_cont_inter_values, unsigned &counter_for_recv_unsigneds, Vector< unsigned > &recv_unsigneds, unsigned &counter_for_recv_doubles, Vector< double > &recv_doubles)
Helper functiono to add external halo node that is not a master.
bool Doc_stats
Boolean to indicate whether to output basic info during setup_multi_domain_interaction() routines.
bool Doc_timings
Boolean to indicate whether to doc timings or not.
Vector< int > Proc_id_plus_one_of_external_element
Proc_id_plus_one_of_external_element[i] contains the processor id (plus one) of the processor on whic...
Definition: multi_domain.cc:91
Vector< double > Flat_packed_located_coordinates
Vector of flat-packed local coordinates for zeta tuples that have been located.
void send_and_receive_located_info(int &iproc, Mesh *const &external_mesh_pt, Problem *problem_pt)
Send location information from current process; Received location information from (current process +...
unsigned Counter_for_flat_packed_unsigneds
Counter used when processing vector of flat-packed unsigneds – this is really "private" data,...
Vector< Vector< unsigned > > External_element_located
Lookup scheme for whether a local element's integration point has had an external element assigned to...
Definition: multi_domain.cc:68
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
Vector< double > Flat_packed_doubles
Vector of flat-packed doubles to be communicated with other processors.
void send_and_receive_missing_zetas(Problem *problem_pt)
Send the zeta coordinates from the current process to the next process; receive from the previous pro...
void setup_bulk_elements_adjacent_to_face_mesh(Problem *problem_pt, Vector< unsigned > &boundary_in_bulk_mesh, Mesh *const &bulk_mesh_pt, Vector< Mesh * > &face_mesh_pt, const unsigned &interaction=0)
Identify the FaceElements (stored in the mesh pointed to by face_mesh_pt) that are adjacent to the bu...
void clean_up()
Helper function that clears all the information used during the external storage creation.
void setup_multi_domain_interactions(Problem *problem_pt, Mesh *const &first_mesh_pt, Mesh *const &second_mesh_pt, const unsigned &first_interaction=0, const unsigned &second_interaction=0)
Set up the two-way multi-domain interactions for the problem pointed to by problem_pt....
Vector< unsigned > Flat_packed_unsigneds
Vector of flat-packed unsigneds to be communicated with other processors – this is really "private" d...
void add_external_halo_node_to_storage(Node *&new_nod_pt, Mesh *const &external_mesh_pt, unsigned &loc_p, unsigned &node_index, FiniteElement *const &new_el_pt, int &n_cont_inter_values, Problem *problem_pt)
Helper function to add external halo nodes, including any masters, based on information received from...
void get_dim_helper(Problem *problem_pt, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt)
Helper function that computes the dimension of the elements within each of the specified meshes (and ...
std::ofstream Doc_boundary_coordinate_file
Output file to document the boundary coordinate along the mesh boundary of the bulk mesh during call ...
Definition: multi_domain.cc:47
Vector< double > Flat_packed_zetas_not_found_locally
Vector of flat-packed zeta coordinates for which the external element could not be found during curre...
Definition: multi_domain.cc:74
void recursively_add_masters_of_external_halo_node_to_storage(Node *&new_nod_pt, Mesh *const &external_mesh_pt, unsigned &loc_p, unsigned &node_index, FiniteElement *const &new_el_pt, int &n_cont_inter_values, Problem *problem_pt)
Recursively add masters of external halo nodes (and their masters, etc) based on information received...
void locate_zeta_for_local_coordinates(const Vector< Mesh * > &mesh_pt, Mesh *const &external_mesh_pt, Vector< MeshAsGeomObject * > &mesh_geom_obj_pt, const unsigned &interaction_index)
locate zeta for current set of "local" coordinates vector-based version
void construct_new_external_halo_master_node_helper(Node *&new_master_nod_pt, Node *&nod_pt, unsigned &loc_p, Mesh *const &external_mesh_pt, Problem *problem_pt)
Helper function which constructs a new external halo master node with the information sent from the h...
void create_external_halo_elements(int &iproc, const Vector< Mesh * > &mesh_pt, Mesh *const &external_mesh_pt, Problem *problem_pt, const unsigned &interaction_index)
Create external (halo) elements on the loop process based on the information received from each locat...
void add_external_halo_master_node_helper(Node *&new_master_nod_pt, Node *&new_nod_pt, Mesh *const &external_mesh_pt, unsigned &loc_p, int &n_cont_inter_values, Problem *problem_pt)
Helper function to add external halo node that is a master.
bool Use_bulk_element_as_external
Boolean to indicate when to use the bulk element as the external element. Defaults to false,...
void setup_multi_domain_interaction(Problem *problem_pt, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, const unsigned &interaction_index=0)
Function to set up the one-way multi-domain interaction for problems where the meshes pointed to by m...
void aux_setup_multi_domain_interaction(Problem *problem_pt, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, const unsigned &interaction_index, Mesh *const &external_face_mesh_pt=0)
Auxiliary helper function.
Vector< unsigned > Located_element_status
Vector to indicate (to another processor) whether a located element (that will have to represented as...
Definition: multi_domain.cc:98
unsigned Counter_for_flat_packed_doubles
Counter used when processing vector of flat-packed doubles – this is really "private" data,...
void locate_zeta_for_missing_coordinates(int &iproc, Mesh *const &external_mesh_pt, Problem *problem_pt, Vector< MeshAsGeomObject * > &mesh_geom_obj_pt)
Locate zeta for current set of missing coordinates; vector-based version.
bool Accept_failed_locate_zeta_in_setup_multi_domain_interaction
Boolean to indicate that failure in setup multi domain functions is acceptable; defaults to false....
Definition: multi_domain.cc:56
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...