problem.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
27#ifdef OOMPH_HAS_MPI
28#include "mpi.h"
29#endif
30
31#include <list>
32#include <algorithm>
33#include <string>
34
35#include "oomph_utilities.h"
36#include "problem.h"
37#include "timesteppers.h"
40#include "refineable_mesh.h"
41#include "triangle_mesh.h"
42#include "linear_solver.h"
43#include "eigen_solver.h"
44#include "assembly_handler.h"
45#include "dg_elements.h"
46#include "partitioning.h"
47#include "spines.h"
48
49// Include to fill in additional_setup_shared_node_scheme() function
51
52
53namespace oomph
54{
55 /// ///////////////////////////////////////////////////////////////
56 // Non-inline functions for the problem class
57 /// ///////////////////////////////////////////////////////////////
58
59 //=================================================================
60 /// The continuation timestepper object
61 //=================================================================
62 ContinuationStorageScheme Problem::Continuation_time_stepper;
63
64 //================================================================
65 /// Constructor: Allocate space for one time stepper
66 /// and set all pointers to NULL and set defaults for all
67 /// parameters.
68 //===============================================================
70 : Mesh_pt(0),
71 Time_pt(0),
72 Explicit_time_stepper_pt(0),
73 Saved_dof_pt(0),
74 Default_set_initial_condition_called(false),
75 Use_globally_convergent_newton_method(false),
76 Empty_actions_before_read_unstructured_meshes_has_been_called(false),
77 Empty_actions_after_read_unstructured_meshes_has_been_called(false),
78 Store_local_dof_pt_in_elements(false),
79 Calculate_hessian_products_analytic(false),
80#ifdef OOMPH_HAS_MPI
81 Doc_imbalance_in_parallel_assembly(false),
82 Use_default_partition_in_load_balance(false),
83 Must_recompute_load_balance_for_assembly(true),
84 Halo_scheme_pt(0),
85#endif
86 Relaxation_factor(1.0),
87 Newton_solver_tolerance(1.0e-8),
89 Nnewton_iter_taken(0),
90 Max_residuals(10.0),
91 Time_adaptive_newton_crash_on_solve_fail(false),
92 Jacobian_reuse_is_enabled(false),
93 Jacobian_has_been_computed(false),
94 Problem_is_nonlinear(true),
95 Pause_at_end_of_sparse_assembly(false),
96 Doc_time_in_distribute(false),
97 Sparse_assembly_method(Perform_assembly_using_vectors_of_pairs),
98 Sparse_assemble_with_arrays_initial_allocation(400),
99 Sparse_assemble_with_arrays_allocation_increment(150),
100 Numerical_zero_for_sparse_assembly(0.0),
101 FD_step_used_in_get_hessian_vector_products(1.0e-8),
102 Mass_matrix_reuse_is_enabled(false),
103 Mass_matrix_has_been_computed(false),
104 Discontinuous_element_formulation(false),
105 Minimum_dt(1.0e-12),
106 Maximum_dt(1.0e12),
107 DTSF_max_increase(4.0),
108 DTSF_min_decrease(0.8),
109 Minimum_dt_but_still_proceed(-1.0),
110 Scale_arc_length(true),
111 Desired_proportion_of_arc_length(0.5),
112 Theta_squared(1.0),
113 Sign_of_jacobian(0),
114 Continuation_direction(1.0),
115 Parameter_derivative(1.0),
116 Parameter_current(0.0),
117 Use_continuation_timestepper(false),
118 Dof_derivative_offset(1),
119 Dof_current_offset(2),
120 Ds_current(0.0),
121 Desired_newton_iterations_ds(5),
122 Minimum_ds(1.0e-10),
123 Bifurcation_detection(false),
124 Bisect_to_find_bifurcation(false),
125 First_jacobian_sign_change(false),
126 Arc_length_step_taken(false),
127 Use_finite_differences_for_continuation_derivatives(false),
128#ifdef OOMPH_HAS_MPI
129 Dist_problem_matrix_distribution(Uniform_matrix_distribution),
130 Parallel_sparse_assemble_previous_allocation(0),
131 Problem_has_been_distributed(false),
132 Bypass_increase_in_dof_check_during_pruning(false),
133 Max_permitted_error_for_halo_check(1.0e-14),
134#endif
135 Shut_up_in_newton_solve(false),
136 Always_take_one_newton_step(false),
137 Timestep_reduction_factor_after_nonconvergence(0.5),
138 Keep_temporal_error_below_tolerance(true)
139 {
141
142 /// Setup terminate helper
144
145 // By default no submeshes:
146 Sub_mesh_pt.resize(0);
147 // No timesteppers
148 Time_stepper_pt.resize(0);
149
150 // Set the linear solvers, eigensolver and assembly handler
153
155
157
158 // setup the communicator
159#ifdef OOMPH_HAS_MPI
161 {
163 }
164 else
165 {
167 }
168#else
170#endif
171
172 // just create an empty linear algebra distribution for the
173 // DOFs
174 // this is setup when assign_eqn_numbers(...) is called.
176 }
177
178 //================================================================
179 /// Destructor to clean up memory
180 //================================================================
182 {
183 // Delete the memory assigned for the global time
184 // (it's created on the fly in Problem::add_time_stepper_pt()
185 // so we are entitled to delete it.
186 if (Time_pt != 0)
187 {
188 delete Time_pt;
189 Time_pt = 0;
190 }
191
192 // We're not using the default linear solver,
193 // somebody else must have built it, so that person
194 // must be in charge of killing it.
195 // We can safely delete the defaults, however
197
200 delete Communicator_pt;
201 delete Dof_distribution_pt;
202
203 // Delete any copies of the problem that have been created for
204 // use in adaptive bifurcation tracking.
205 // ALH: This will eventually go
206 unsigned n_copies = Copy_of_problem_pt.size();
207 for (unsigned c = 0; c < n_copies; c++)
208 {
209 delete Copy_of_problem_pt[c];
210 }
211
212 // if this problem has sub meshes then we must delete the Mesh_pt
213 if (Sub_mesh_pt.size() != 0)
214 {
216 delete Mesh_pt;
217 }
218
219 // Since we called the TerminateHelper setup function in the constructor,
220 // we need to delete anything that was dynamically allocated (as it's
221 // just a namespace and so doesn't have it's own destructor) in the function
223 }
224
225 //=================================================================
226 /// Setup the count vector that records how many elements contribute
227 /// to each degree of freedom. Returns the total number of elements
228 /// in the problem
229 //=================================================================
231 {
232 // Now set the element counter to have the current Dof distribution
234 // We need to use the halo scheme (assuming it has been setup)
235#ifdef OOMPH_HAS_MPI
237#endif
238
239 // Loop over the elements and count the entries
240 // and number of (non-halo) elements
241 const unsigned n_element = this->mesh_pt()->nelement();
242 unsigned n_non_halo_element_local = 0;
243 for (unsigned e = 0; e < n_element; e++)
244 {
245 GeneralisedElement* elem_pt = this->mesh_pt()->element_pt(e);
246#ifdef OOMPH_HAS_MPI
247 // Ignore halo elements
248 if (!elem_pt->is_halo())
249 {
250#endif
251 // Increment the number of non halo elements
252 ++n_non_halo_element_local;
253 // Now count the number of times the element contributes to a value
254 // using the current assembly handler
255 unsigned n_var = this->Assembly_handler_pt->ndof(elem_pt);
256 for (unsigned n = 0; n < n_var; n++)
257 {
259 this->Assembly_handler_pt->eqn_number(elem_pt, n));
260 }
261#ifdef OOMPH_HAS_MPI
262 }
263#endif
264 }
265
266 // Storage for the total number of elements
267 unsigned Nelement = 0;
268
269 // Add together all the counts if we are in parallel
270#ifdef OOMPH_HAS_MPI
272
273 // If distributed, find the total number of elements in the problem
275 {
276 // Need to gather the total number of non halo elements
277 MPI_Allreduce(&n_non_halo_element_local,
278 &Nelement,
279 1,
280 MPI_UNSIGNED,
281 MPI_SUM,
282 this->communicator_pt()->mpi_comm());
283 }
284 // Otherwise the total number is the same on each processor
285 else
286#endif
287 {
288 Nelement = n_non_halo_element_local;
289 }
290
291 return Nelement;
292 }
293
294
295#ifdef OOMPH_HAS_MPI
296
297 //==================================================================
298 /// Setup the halo scheme for the degrees of freedom
299 //==================================================================
301 {
302 // Find the number of elements stored on this processor
303 const unsigned n_element = this->mesh_pt()->nelement();
304
305 // Work out the all global equations to which this processor
306 // contributes
307 Vector<unsigned> my_eqns;
308 this->get_my_eqns(this->Assembly_handler_pt, 0, n_element - 1, my_eqns);
309
310 // Build the halo scheme, based on the equations to which this
311 // processor contributes
314
315 // Find pointers to all the halo dofs
316 // There may be more of these than required by my_eqns
317 //(but should not be less)
318 std::map<unsigned, double*> halo_data_pt;
319 this->get_all_halo_data(halo_data_pt);
320
321 // Now setup the Halo_dofs
322 Halo_scheme_pt->setup_halo_dofs(halo_data_pt, this->Halo_dof_pt);
323 }
324
325 //==================================================================
326 /// Distribute the problem without doc; report stats if required.
327 /// Returns actual partitioning used, e.g. for restart.
328 //==================================================================
329 Vector<unsigned> Problem::distribute(const bool& report_stats)
330 {
331 // Set dummy doc paramemters
332 DocInfo doc_info;
333 doc_info.disable_doc();
334
335 // Set the sizes of the input and output vectors
336 unsigned n_element = mesh_pt()->nelement();
337 Vector<unsigned> element_partition(n_element, 0);
338
339 // Distribute and return partitioning
340 return distribute(element_partition, doc_info, report_stats);
341 }
342
343 //==================================================================
344 /// Distribute the problem according to specified partition.
345 /// If all entries in partitioning vector are zero we use METIS
346 /// to do the partitioning after all.
347 /// Returns actual partitioning used, e.g. for restart.
348 //==================================================================
350 const Vector<unsigned>& element_partition, const bool& report_stats)
351 {
352#ifdef PARANOID
353 bool has_non_zero_entry = false;
354 unsigned n = element_partition.size();
355 for (unsigned i = 0; i < n; i++)
356 {
357 if (element_partition[i] != 0)
358 {
359 has_non_zero_entry = true;
360 break;
361 }
362 }
363 if (!has_non_zero_entry)
364 {
365 std::ostringstream warn_message;
366 warn_message << "WARNING: All entries in specified partitioning vector \n"
367 << " are zero -- will ignore this and use METIS\n"
368 << " to perform the partitioning\n";
370 warn_message.str(), "Problem::distribute()", OOMPH_EXCEPTION_LOCATION);
371 }
372#endif
373 // Set dummy doc paramemters
374 DocInfo doc_info;
375 doc_info.disable_doc();
376
377 // Distribute and return partitioning
378 return distribute(element_partition, doc_info, report_stats);
379 }
380
381 //==================================================================
382 /// Distribute the problem and doc to specified DocInfo.
383 /// Returns actual partitioning used, e.g. for restart.
384 //==================================================================
386 const bool& report_stats)
387 {
388 // Set the sizes of the input and output vectors
389 unsigned n_element = mesh_pt()->nelement();
390
391 // Dummy input vector
392 Vector<unsigned> element_partition(n_element, 0);
393
394 // Distribute and return partitioning
395 return distribute(element_partition, doc_info, report_stats);
396 }
397
398 //==================================================================
399 /// Distribute the problem according to specified partition.
400 /// (If all entries in partitioning vector are zero we use METIS
401 /// to do the partitioning after all) and doc.
402 /// Returns actual partitioning used, e.g. for restart.
403 //==================================================================
405 const Vector<unsigned>& element_partition,
406 DocInfo& doc_info,
407 const bool& report_stats)
408 {
409 // Storage for number of processors and number of elements in global mesh
410 int n_proc = this->communicator_pt()->nproc();
411 int my_rank = this->communicator_pt()->my_rank();
412 int n_element = mesh_pt()->nelement();
413
414 // Vector to be returned
415 Vector<unsigned> return_element_domain;
416
417 // Buffer extreme cases
418 if (n_proc == 1) // single-process job - don't do anything
419 {
420 if (report_stats)
421 {
422 std::ostringstream warn_message;
423 warn_message << "WARNING: You've tried to distribute a problem over\n"
424 << "only one processor: this would make METIS crash.\n"
425 << "Ignoring your request for distribution.\n";
426 OomphLibWarning(warn_message.str(),
427 "Problem::distribute()",
428 OOMPH_EXCEPTION_LOCATION);
429 }
430 }
431 else if (n_proc > n_element) // more processors than elements
432 {
433 // Throw an error
434 std::ostringstream error_stream;
435 error_stream << "You have tried to distribute a problem\n"
436 << "but there are less elements than processors.\n"
437 << "Please re-run with more elements!\n"
438 << "Please also ensure that actions_before_distribute().\n"
439 << "and actions_after_distribute() are correctly set up.\n"
440 << std::endl;
441 throw OomphLibError(
442 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
443 }
444 else
445 {
446 // We only distribute uniformly-refined meshes; buffer the case where
447 // either mesh is not uniformly refined
448 bool a_mesh_is_not_uniformly_refined = false;
449 unsigned n_mesh = nsub_mesh();
450 if (n_mesh == 0)
451 {
452 // Check refinement levels
453 if (TreeBasedRefineableMeshBase* mmesh_pt =
454 dynamic_cast<TreeBasedRefineableMeshBase*>(mesh_pt(0)))
455 {
456 unsigned min_ref_level = 0;
457 unsigned max_ref_level = 0;
458 mmesh_pt->get_refinement_levels(min_ref_level, max_ref_level);
459 // If they are not the same
460 if (max_ref_level != min_ref_level)
461 {
462 a_mesh_is_not_uniformly_refined = true;
463 }
464 }
465 }
466 else
467 {
468 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
469 {
470 // Check refinement levels for each mesh individually
471 // (one mesh is allowed to be "more uniformly refined" than another)
472 if (TreeBasedRefineableMeshBase* mmesh_pt =
473 dynamic_cast<TreeBasedRefineableMeshBase*>(mesh_pt(i_mesh)))
474 {
475 unsigned min_ref_level = 0;
476 unsigned max_ref_level = 0;
477 mmesh_pt->get_refinement_levels(min_ref_level, max_ref_level);
478 // If they are not the same
479 if (max_ref_level != min_ref_level)
480 {
481 a_mesh_is_not_uniformly_refined = true;
482 }
483 }
484 }
485 }
486
487 // If any mesh is not uniformly refined
488 if (a_mesh_is_not_uniformly_refined)
489 {
490 // Again it may make more sense to throw an error here as the user
491 // will probably not be running a problem that is small enough to
492 // fit the whole of on each processor
493 std::ostringstream error_stream;
494 error_stream << "You have tried to distribute a problem\n"
495 << "but at least one of your meshes is no longer\n"
496 << "uniformly refined. In order to preserve the Tree\n"
497 << "and TreeForest structure, Problem::distribute() can\n"
498 << "only be called while meshes are uniformly refined.\n"
499 << std::endl;
500 throw OomphLibError(
501 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
502 }
503 else
504 {
505 // Is there any global data? If so, distributing the problem won't work
506 if (nglobal_data() > 0)
507 {
508 std::ostringstream error_stream;
509 error_stream << "You have tried to distribute a problem\n"
510 << "and there is some global data.\n"
511 << "This is not likely to work...\n"
512 << std::endl;
513 throw OomphLibError(error_stream.str(),
514 OOMPH_CURRENT_FUNCTION,
515 OOMPH_EXCEPTION_LOCATION);
516 }
517
518 double t_start = 0;
520 {
521 t_start = TimingHelpers::timer();
522 }
523
524
525#ifdef PARANOID
526 unsigned old_ndof = ndof();
527#endif
528
529 // Need to partition the global mesh before distributing
530 Mesh* global_mesh_pt = mesh_pt();
531
532 // Vector listing the affiliation of each element
533 unsigned nelem = global_mesh_pt->nelement();
534 Vector<unsigned> element_domain(nelem);
535
536 // Number of elements that I'm in charge of, based on any
537 // incoming partitioning
538 unsigned n_my_elements = 0;
539
540 // Have we used the pre-set partitioning
541 bool used_preset_partitioning = false;
542
543 // Partition the mesh, unless the partition has already been passed in
544 // If it hasn't then the sum of all the entries of the vector should be
545 // 0
546 unsigned sum_element_partition = 0;
547 unsigned n_part = element_partition.size();
548 for (unsigned e = 0; e < n_part; e++)
549 {
550 // ... another one for me.
551 if (int(element_partition[e]) == my_rank) n_my_elements++;
552
553 sum_element_partition += element_partition[e];
554 }
555 if (sum_element_partition == 0)
556 {
557 oomph_info << "INFO: using METIS to partition elements" << std::endl;
558 partition_global_mesh(global_mesh_pt, doc_info, element_domain);
559 used_preset_partitioning = false;
560 }
561 else
562 {
563 oomph_info << "INFO: using pre-set partition of elements"
564 << std::endl;
565 used_preset_partitioning = true;
566 element_domain = element_partition;
567 }
568
569 // Set the GLOBAL Mesh as being distributed
570 global_mesh_pt->set_communicator_pt(this->communicator_pt());
571
572 double t_end = 0.0;
574 {
575 t_end = TimingHelpers::timer();
576 oomph_info << "Time for partitioning of global mesh: "
577 << t_end - t_start << std::endl;
578 t_start = TimingHelpers::timer();
579 }
580
581 // Store how many elements we had in the various sub-meshes
582 // before actions_before_distribute() (which may empty some of
583 // them).
584 Vector<unsigned> n_element_in_old_submesh(n_mesh);
585 if (n_mesh != 0)
586 {
587 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
588 {
589 unsigned nsub_elem = mesh_pt(i_mesh)->nelement();
590 n_element_in_old_submesh[i_mesh] = nsub_elem;
591 }
592 }
593
594 // Partitioning complete; call actions before distribute
596
598 {
599 t_end = TimingHelpers::timer();
600 oomph_info << "Time for actions before distribute: "
601 << t_end - t_start << std::endl;
602 }
603
604 // This next bit is cheap -- omit timing
605 // t_start = TimingHelpers::timer();
606
607 // Number of submeshes (NB: some may have been deleted in
608 // actions_after_distribute())
609 n_mesh = nsub_mesh();
610
611
612 // Prepare vector of vectors for submesh element domains
613 Vector<Vector<unsigned>> submesh_element_domain(n_mesh);
614
615 // The submeshes need to know their own element domains.
616 // Also if any meshes have been emptied we ignore their
617 // partitioning in the vector that we return from here
618 return_element_domain.reserve(element_domain.size());
619 if (n_mesh != 0)
620 {
621 unsigned count = 0;
622 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
623 {
624 unsigned nsub_elem = mesh_pt(i_mesh)->nelement();
625 submesh_element_domain[i_mesh].resize(nsub_elem);
626 unsigned nsub_elem_old = n_element_in_old_submesh[i_mesh];
627 for (unsigned e = 0; e < nsub_elem_old; e++)
628 {
629 if (nsub_elem_old == nsub_elem)
630 {
631 submesh_element_domain[i_mesh][e] = element_domain[count];
632 return_element_domain.push_back(element_domain[count]);
633 }
634 // return_element_domain.push_back(element_domain[count]);
635 count++;
636 }
637 }
638 }
639 else
640 {
641 return_element_domain = element_domain;
642 }
643
645 {
646 t_start = TimingHelpers::timer();
647 }
648
649 // Setup the map between "root" element and number in global mesh
650 // (currently used in the load_balance() routines)
651
652 // This map is only established for structured meshes, then we
653 // need to check here the type of mesh
654 if (n_mesh == 0)
655 {
656 // Check if the only one mesh is an structured mesh
657 bool structured_mesh = true;
658 TriangleMeshBase* tri_mesh_pt =
659 dynamic_cast<TriangleMeshBase*>(mesh_pt(0));
660 if (tri_mesh_pt != 0)
661 {
662 structured_mesh = false;
663 } // if (tri_mesh_pt != 0)
664 if (structured_mesh)
665 {
666 const unsigned n_ele = global_mesh_pt->nelement();
667 Base_mesh_element_pt.resize(n_ele);
669 for (unsigned e = 0; e < n_ele; e++)
670 {
671 GeneralisedElement* el_pt = global_mesh_pt->element_pt(e);
673 Base_mesh_element_pt[e] = el_pt;
674 } // for (e<n_ele)
675 } // A TreeBaseMesh mesh
676 } // if (n_mesh==0)
677 else
678 {
679 // If we have submeshes then we only add those elements that
680 // belong to structured meshes, but first compute the number
681 // of total elements in the structured meshes
682 unsigned nglobal_element = 0;
683 // Store which submeshes are structured
684 std::vector<bool> is_structured_mesh(n_mesh);
685 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
686 {
687 TriangleMeshBase* tri_mesh_pt =
688 dynamic_cast<TriangleMeshBase*>(mesh_pt(i_mesh));
689 if (tri_mesh_pt != 0)
690 {
691 // Set the flags to indicate this is not an structured
692 // mesh
693 is_structured_mesh[i_mesh] = false;
694 } // if (tri_mesh_pt != 0)
695 else
696 {
697 // Set the flags to indicate this is an structured
698 // mesh
699 is_structured_mesh[i_mesh] = true;
700 } // else if (tri_mesh_pt != 0)
701 // Check if mesh is an structured mesh
702 if (is_structured_mesh[i_mesh])
703 {
704 nglobal_element += mesh_pt(i_mesh)->nelement();
705 } // A TreeBaseMesh mesh
706 } // for (i_mesh<n_mesh)
707
708 // Once computed the number of elements, then resize the
709 // structure
710 Base_mesh_element_pt.resize(nglobal_element);
712 unsigned counter = 0;
713 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
714 {
715 // Check if mesh is an structured mesh
716 if (is_structured_mesh[i_mesh])
717 {
718 const unsigned n_ele = mesh_pt(i_mesh)->nelement();
719 for (unsigned e = 0; e < n_ele; e++)
720 {
721 GeneralisedElement* el_pt = mesh_pt(i_mesh)->element_pt(e);
722 Base_mesh_element_number_plus_one[el_pt] = counter + 1;
723 Base_mesh_element_pt[counter] = el_pt;
724 // Inrease the global element number
725 counter++;
726 } // for (e<n_ele)
727 } // An structured mesh
728 } // for (i_mesh<n_mesh)
729
730#ifdef PARANOID
731 if (counter != nglobal_element)
732 {
733 std::ostringstream error_stream;
734 error_stream
735 << "The number of global elements (" << nglobal_element
736 << ") is not the sameas the number of\nadded elements ("
737 << counter << ") to the Base_mesh_element_pt data "
738 << "structure!!!\n\n";
739 throw OomphLibError(error_stream.str(),
740 "Problem::distribute()",
741 OOMPH_EXCEPTION_LOCATION);
742 } // if (counter != nglobal_element)
743#endif // #ifdef PARANOID
744
745 } // else if (n_mesh==0)
746
747 // Wipe everything if a pre-determined partitioning
748 // didn't specify ANY elements for this processor
749 // (typically happens during restarts with larger number
750 // of processors -- in this case we really want an empty
751 // processor rather than one with any "kept" halo elements)
752 bool overrule_keep_as_halo_element_status = false;
753 if ((n_my_elements == 0) && (used_preset_partitioning))
754 {
755 oomph_info << "INFO: We're over-ruling the \"keep as halo element\"\n"
756 << " status because the preset partitioning\n"
757 << " didn't place ANY elements on this processor,\n"
758 << " probably because of a restart on a larger \n"
759 << " number of processors\n";
760 overrule_keep_as_halo_element_status = true;
761 }
762
763
764 // Distribute the (sub)meshes (i.e. sort out their halo lookup schemes)
765 Vector<GeneralisedElement*> deleted_element_pt;
766 if (n_mesh == 0)
767 {
768 global_mesh_pt->distribute(this->communicator_pt(),
769 element_domain,
770 deleted_element_pt,
771 doc_info,
772 report_stats,
773 overrule_keep_as_halo_element_status);
774 }
775 else // There are submeshes, "distribute" each one separately
776 {
777 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
778 {
779 if (report_stats)
780 {
781 oomph_info << "Distributing submesh " << i_mesh << std::endl
782 << "--------------------" << std::endl;
783 }
784 // Set the doc_info number to reflect the submesh
785 doc_info.number() = i_mesh;
786 mesh_pt(i_mesh)->distribute(this->communicator_pt(),
787 submesh_element_domain[i_mesh],
788 deleted_element_pt,
789 doc_info,
790 report_stats,
791 overrule_keep_as_halo_element_status);
792 }
793 // Rebuild the global mesh
795 }
796
797 // Null out information associated with deleted elements
798 unsigned n_del = deleted_element_pt.size();
799 for (unsigned e = 0; e < n_del; e++)
800 {
801 GeneralisedElement* el_pt = deleted_element_pt[e];
802 unsigned old_el_number = Base_mesh_element_number_plus_one[el_pt] - 1;
804 Base_mesh_element_pt[old_el_number] = 0;
805 }
806
808 {
809 t_end = TimingHelpers::timer();
810 oomph_info << "Time for mesh-level distribution: " << t_end - t_start
811 << std::endl;
812 t_start = TimingHelpers::timer();
813 }
814
815 // Now the problem has been distributed
817
818 // Call actions after distribute
820
822 {
823 t_end = TimingHelpers::timer();
824 oomph_info << "Time for actions after distribute: " << t_end - t_start
825 << std::endl;
826 t_start = TimingHelpers::timer();
827 }
828
829 // Re-assign the equation numbers (incl synchronisation if reqd)
830 unsigned n_dof = assign_eqn_numbers();
831 oomph_info << "Number of equations: " << n_dof << std::endl;
832
834 {
835 t_end = TimingHelpers::timer();
836 oomph_info << "Time for re-assigning eqn numbers (in distribute): "
837 << t_end - t_start << std::endl;
838 }
839
840
841#ifdef PARANOID
842 if (n_dof != old_ndof)
843 {
844 std::ostringstream error_stream;
845 error_stream
846 << "Number of dofs in distribute() has changed "
847 << "from " << old_ndof << " to " << n_dof << "\n"
848 << "Check that you've implemented any necessary "
849 "actions_before/after\n"
850 << "distribute functions, e.g. to pin redundant pressure dofs"
851 << " etc.\n";
852 throw OomphLibError(error_stream.str(),
853 OOMPH_CURRENT_FUNCTION,
854 OOMPH_EXCEPTION_LOCATION);
855 }
856#endif
857
858 } // end if to check for uniformly refined mesh(es)
859
860 } // end if to check number of processors vs. number of elements etc.
861
862
863 // Force re-analysis of time spent on assembly each
864 // elemental Jacobian
867
868 // Return the partition vector used in the distribution
869 return return_element_domain;
870 }
871
872 //==================================================================
873 /// Partition the global mesh, return vector specifying the processor
874 /// number for each element. Virtual so that it can be overloaded by
875 /// any user; the default is to use METIS to perform the partitioning
876 /// (with a bit of cleaning up afterwards to sort out "special cases").
877 //==================================================================
879 DocInfo& doc_info,
880 Vector<unsigned>& element_domain,
881 const bool& report_stats)
882 {
883 // Storage for number of processors and current processor
884 int n_proc = this->communicator_pt()->nproc();
885 int rank = this->communicator_pt()->my_rank();
886
887 std::ostringstream filename;
888 std::ofstream some_file;
889
890 // Doc the original mesh on proc 0
891 //--------------------------------
892 if (doc_info.is_doc_enabled())
893 {
894 if (rank == 0)
895 {
896 filename << doc_info.directory() << "/complete_mesh"
897 << doc_info.number() << ".dat";
898 global_mesh_pt->output(filename.str().c_str(), 5);
899 }
900 }
901
902 // Partition the mesh
903 //-------------------
904 // METIS Objective (0: minimise edge cut; 1: minimise total comm volume)
905 unsigned objective = 0;
906
907 // Do the partitioning
908 unsigned nelem = 0;
909 if (this->communicator_pt()->my_rank() == 0)
910 {
911 METIS::partition_mesh(this, n_proc, objective, element_domain);
912 nelem = element_domain.size();
913 }
914 MPI_Bcast(&nelem, 1, MPI_UNSIGNED, 0, this->communicator_pt()->mpi_comm());
915 element_domain.resize(nelem);
916 MPI_Bcast(&element_domain[0],
917 nelem,
918 MPI_UNSIGNED,
919 0,
920 this->communicator_pt()->mpi_comm());
921
922 // On very coarse meshes with larger numbers of processors, METIS
923 // occasionally returns an element_domain Vector for which a particular
924 // processor has no elements affiliated to it; the following fixes this
925
926 // Convert element_domain to integer storage
927 Vector<int> int_element_domain(nelem);
928 for (unsigned e = 0; e < nelem; e++)
929 {
930 int_element_domain[e] = element_domain[e];
931 }
932
933 // Global storage for number of elements on each process
934 int my_number_of_elements = 0;
935 Vector<int> number_of_elements(n_proc, 0);
936
937 for (unsigned e = 0; e < nelem; e++)
938 {
939 if (int_element_domain[e] == rank)
940 {
941 my_number_of_elements++;
942 }
943 }
944
945 // Communicate the correct value for each single process into
946 // the global storage vector
947 MPI_Allgather(&my_number_of_elements,
948 1,
949 MPI_INT,
950 &number_of_elements[0],
951 1,
952 MPI_INT,
953 this->communicator_pt()->mpi_comm());
954
955 // If a process has no elements then switch an element with the
956 // process with the largest number of elements, assuming
957 // that it still has enough elements left to share
958 int max_number_of_elements = 0;
959 int process_with_max_elements = 0;
960 for (int d = 0; d < n_proc; d++)
961 {
962 if (number_of_elements[d] == 0)
963 {
964 // Find the process with maximum number of elements
965 if (max_number_of_elements <= 1)
966 {
967 for (int dd = 0; dd < n_proc; dd++)
968 {
969 if (number_of_elements[dd] > max_number_of_elements)
970 {
971 max_number_of_elements = number_of_elements[dd];
972 process_with_max_elements = dd;
973 }
974 }
975 }
976
977 // Check that this number of elements is okay for sharing
978 if (max_number_of_elements <= 1)
979 {
980 // Throw an error if elements can't be shared
981 std::ostringstream error_stream;
982 error_stream << "No process has more than 1 element, and\n"
983 << "at least one process has no elements!\n"
984 << "Suggest rerunning with more refinement.\n"
985 << std::endl;
986 throw OomphLibError(error_stream.str(),
987 OOMPH_CURRENT_FUNCTION,
988 OOMPH_EXCEPTION_LOCATION);
989 }
990
991 // Loop over the element domain vector and switch
992 // one value for process "process_with_max_elements" with d
993 for (unsigned e = 0; e < nelem; e++)
994 {
995 if (int_element_domain[e] == process_with_max_elements)
996 {
997 int_element_domain[e] = d;
998 // Change the numbers associated with these processes
999 number_of_elements[d]++;
1000 number_of_elements[process_with_max_elements]--;
1001 // Reduce the number of elements available on "max" process
1002 max_number_of_elements--;
1003 // Inform the user that a switch has taken place
1004 if (report_stats)
1005 {
1006 oomph_info << "INFO: Switched element domain at position " << e
1007 << std::endl
1008 << "from process " << process_with_max_elements
1009 << " to process " << d << std::endl
1010 << "which was given no elements by METIS partition"
1011 << std::endl;
1012 }
1013 // Only need to do this once for this element loop, otherwise
1014 // this will take all the elements from "max" process and put them
1015 // in process d, thus leaving essentially the same problem!
1016 break;
1017 }
1018 }
1019 }
1020 }
1021
1022 // Reassign new values to the element_domain vector
1023 for (unsigned e = 0; e < nelem; e++)
1024 {
1025 element_domain[e] = int_element_domain[e];
1026 }
1027
1028 unsigned count_elements = 0;
1029 for (unsigned e = 0; e < nelem; e++)
1030 {
1031 if (int(element_domain[e]) == rank)
1032 {
1033 count_elements++;
1034 }
1035 }
1036
1037 if (report_stats)
1038 {
1039 oomph_info << "I have " << count_elements
1040 << " elements from this partition" << std::endl
1041 << std::endl;
1042 }
1043 }
1044
1045 //==================================================================
1046 /// (Irreversibly) prune halo(ed) elements and nodes, usually
1047 /// after another round of refinement, to get rid of
1048 /// excessively wide halo layers. Note that the current
1049 /// mesh will be now regarded as the base mesh and no unrefinement
1050 /// relative to it will be possible once this function
1051 /// has been called.
1052 //==================================================================
1054 const bool& report_stats)
1055 {
1056 // Storage for number of processors and current processor
1057 int n_proc = this->communicator_pt()->nproc();
1058
1059 // Has the problem been distributed yet?
1061 {
1063 << "WARNING: Problem::prune_halo_elements_and_nodes() was called on a "
1064 << "non-distributed Problem!" << std::endl;
1065 oomph_info << "Ignoring your request..." << std::endl;
1066 }
1067 else
1068 {
1069 // There are no halo layers to prune if it's a single-process job
1070 if (n_proc == 1)
1071 {
1073 << "WARNING: You've tried to prune halo layers on a problem\n"
1074 << "with only one processor: this is unnecessary.\n"
1075 << "Ignoring your request." << std::endl
1076 << std::endl;
1077 }
1078 else
1079 {
1080#ifdef PARANOID
1081 unsigned old_ndof = ndof();
1082#endif
1083
1084 double t_start = 0.0;
1086 {
1087 t_start = TimingHelpers::timer();
1088 }
1089
1090 // Call actions before distribute
1092
1093 double t_end = 0.0;
1095 {
1096 t_end = TimingHelpers::timer();
1097 oomph_info << "Time for actions_before_distribute() in "
1098 << "Problem::prune_halo_elements_and_nodes(): "
1099 << t_end - t_start << std::endl;
1100 t_start = TimingHelpers::timer();
1101 }
1102
1103 // Associate all elements with root in current Base mesh
1104 unsigned nel = Base_mesh_element_pt.size();
1105 std::map<GeneralisedElement*, unsigned>
1106 old_base_element_number_plus_one;
1107 std::vector<bool> old_root_is_halo_or_non_existent(nel, true);
1108 for (unsigned e = 0; e < nel; e++)
1109 {
1110 // Get the base element
1112
1113 // Does it exist locally?
1114 if (base_el_pt != 0)
1115 {
1116 // Check if it's a halo element
1117 if (!base_el_pt->is_halo())
1118 {
1119 old_root_is_halo_or_non_existent[e] = false;
1120 }
1121
1122 // Not refineable: It's only the element iself
1123 RefineableElement* ref_el_pt = 0;
1124 ref_el_pt = dynamic_cast<RefineableElement*>(base_el_pt);
1125 if (ref_el_pt == 0)
1126 {
1127 old_base_element_number_plus_one[base_el_pt] = e + 1;
1128 }
1129 // Refineable: Get entire tree of elements
1130 else
1131 {
1132 Vector<Tree*> tree_pt;
1133 ref_el_pt->tree_pt()->stick_all_tree_nodes_into_vector(tree_pt);
1134 unsigned ntree = tree_pt.size();
1135 for (unsigned t = 0; t < ntree; t++)
1136 {
1137 old_base_element_number_plus_one[tree_pt[t]->object_pt()] =
1138 e + 1;
1139 }
1140 }
1141 }
1142 }
1143
1144
1146 {
1147 t_end = TimingHelpers::timer();
1148 oomph_info << "Time for setup old root elements in "
1149 << "Problem::prune_halo_elements_and_nodes(): "
1150 << t_end - t_start << std::endl;
1151 t_start = TimingHelpers::timer();
1152 }
1153
1154
1155 // Now remember the old number of base elements
1156 unsigned nel_base_old = nel;
1157
1158
1159 // Prune the halo elements and nodes of the mesh(es)
1160 Vector<GeneralisedElement*> deleted_element_pt;
1161 unsigned n_mesh = nsub_mesh();
1162 if (n_mesh == 0)
1163 {
1164 // Prune halo elements and nodes for the (single) global mesh
1166 deleted_element_pt, doc_info, report_stats);
1167 }
1168 else
1169 {
1170 // Loop over individual submeshes and prune separately
1171 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
1172 {
1174 deleted_element_pt, doc_info, report_stats);
1175 }
1176
1177 // Rebuild the global mesh
1179 }
1180
1182 {
1183 t_end = TimingHelpers::timer();
1184 oomph_info << "Total time for all mesh-level prunes in "
1185 << "Problem::prune_halo_elements_and_nodes(): "
1186 << t_end - t_start << std::endl;
1187 t_start = TimingHelpers::timer();
1188 }
1189
1190 // Loop over all elements in newly rebuilt mesh (which contains
1191 // all element in "tree order"), find the roots
1192 // (which are either non-refineable elements or refineable elements
1193 // whose tree representations are TreeRoots)
1194 std::map<FiniteElement*, bool> root_el_done;
1195
1196 // Vector storing vectors of pointers to new base elements associated
1197 // with the same old base element
1199 new_base_element_associated_with_old_base_element(nel_base_old);
1200
1201 unsigned n_meshes = n_mesh;
1202 // Change the value for the number of submeshes if there is only
1203 // one mesh so that the loop below works if we have only one
1204 // mesh
1205 if (n_meshes == 0)
1206 {
1207 n_meshes = 1;
1208 }
1209
1210 // Store which submeshes, if there are some are structured
1211 // meshes
1212 std::vector<bool> is_structured_mesh(n_meshes);
1213
1214 // Loop over all elements in the rebuilt mesh, but make sure
1215 // that we are only looping over the structured meshes
1216 nel = 0;
1217 for (unsigned i_mesh = 0; i_mesh < n_meshes; i_mesh++)
1218 {
1219 TriangleMeshBase* tri_mesh_pt =
1220 dynamic_cast<TriangleMeshBase*>(mesh_pt(i_mesh));
1221 if (!(tri_mesh_pt != 0))
1222 {
1223 // Mark the mesh as structured mesh
1224 is_structured_mesh[i_mesh] = true;
1225 // Add the number of elements
1226 nel += mesh_pt(i_mesh)->nelement();
1227 } // if (!(tri_mesh_pt!=0))
1228 else
1229 {
1230 // Mark the mesh as nonstructured mesh
1231 is_structured_mesh[i_mesh] = false;
1232 } // else if (!(tri_mesh_pt!=0))
1233 } // for (i_mesh < n_mesh)
1234
1235 // Go for all the meshes (if there are submeshes)
1236 for (unsigned i_mesh = 0; i_mesh < n_meshes; i_mesh++)
1237 {
1238 // Only work with the elements in the mesh if it is an
1239 // structured mesh
1240 if (is_structured_mesh[i_mesh])
1241 {
1242 // Get the number of elements in the submesh
1243 const unsigned nele_submesh = mesh_pt(i_mesh)->nelement();
1244 for (unsigned e = 0; e < nele_submesh; e++)
1245 {
1246 // Get the element
1247 GeneralisedElement* el_pt = mesh_pt(i_mesh)->element_pt(e);
1248
1249 // Not refineable: It's definitely a new base element
1250 RefineableElement* ref_el_pt = 0;
1251 ref_el_pt = dynamic_cast<RefineableElement*>(el_pt);
1252 if (ref_el_pt == 0)
1253 {
1254 unsigned old_base_el_no =
1255 old_base_element_number_plus_one[el_pt] - 1;
1256 new_base_element_associated_with_old_base_element
1257 [old_base_el_no]
1258 .push_back(el_pt);
1259 }
1260 // Refineable
1261 else
1262 {
1263 // Is it a tree root (after pruning)? In that case it's
1264 // a new base element
1265 if (dynamic_cast<TreeRoot*>(ref_el_pt->tree_pt()))
1266 {
1267 unsigned old_base_el_no =
1268 old_base_element_number_plus_one[el_pt] - 1;
1269 new_base_element_associated_with_old_base_element
1270 [old_base_el_no]
1271 .push_back(el_pt);
1272 }
1273 else
1274 {
1275 // Get associated root element
1276 FiniteElement* root_el_pt =
1277 ref_el_pt->tree_pt()->root_pt()->object_pt();
1278
1279 if (!root_el_done[root_el_pt])
1280 {
1281 root_el_done[root_el_pt] = true;
1282 unsigned old_base_el_no =
1283 old_base_element_number_plus_one[el_pt] - 1;
1284 new_base_element_associated_with_old_base_element
1285 [old_base_el_no]
1286 .push_back(root_el_pt);
1287 }
1288 }
1289 }
1290 } // for (e < nele_submesh)
1291 } // if (is_structured_mesh[i_mesh])
1292 } // for (i_mesh < n_mesh)
1293
1294 // Create a vector that stores how many new root/base elements
1295 // got spawned from each old root/base element in the global mesh
1296 Vector<unsigned> local_n_new_root(nel_base_old);
1297#ifdef PARANOID
1298 Vector<unsigned> n_new_root_back(nel_base_old);
1299#endif
1300 for (unsigned e = 0; e < nel_base_old; e++)
1301 {
1302 local_n_new_root[e] =
1303 new_base_element_associated_with_old_base_element[e].size();
1304
1305#ifdef PARANOID
1306 // Backup so we can check that halo data was consistent
1307 n_new_root_back[e] = local_n_new_root[e];
1308#endif
1309 }
1310
1312 {
1313 t_end = TimingHelpers::timer();
1314 oomph_info << "Time for setup of new base elements in "
1315 << "Problem::prune_halo_elements_and_nodes(): "
1316 << t_end - t_start << std::endl;
1317 t_start = TimingHelpers::timer();
1318 }
1319
1320 // Now do reduce operation to get information for all
1321 // old root/base elements -- the pruned (halo!) base elements contain
1322 // fewer associated new roots.
1323 Vector<unsigned> n_new_root(nel_base_old);
1324 MPI_Allreduce(&local_n_new_root[0],
1325 &n_new_root[0],
1326 nel_base_old,
1327 MPI_UNSIGNED,
1328 MPI_MAX,
1329 this->communicator_pt()->mpi_comm());
1330
1331
1333 {
1334 t_end = TimingHelpers::timer();
1335 oomph_info << "Time for allreduce in "
1336 << "Problem::prune_halo_elements_and_nodes(): "
1337 << t_end - t_start << std::endl;
1338 t_start = TimingHelpers::timer();
1339 }
1340
1341 // Find out total number of base elements
1342 unsigned nel_base_new = 0;
1343 for (unsigned e = 0; e < nel_base_old; e++)
1344 {
1345 // Increment
1346 nel_base_new += n_new_root[e];
1347
1348#ifdef PARANOID
1349 // If we already had data for this root previously then
1350 // the data ought to be consistent afterwards (since taking
1351 // the max of consistent numbers shouldn't change things -- this
1352 // deals with halo/haloed elements)
1353 if (!old_root_is_halo_or_non_existent[e])
1354 {
1355 if (n_new_root_back[e] != 0)
1356 {
1357 if (n_new_root_back[e] != n_new_root[e])
1358 {
1359 std::ostringstream error_stream;
1360 error_stream
1361 << "Number of new root elements spawned from old root " << e
1362 << ": " << n_new_root[e] << "\nis not consistent"
1363 << " with previous value: " << n_new_root_back[e]
1364 << std::endl;
1365 throw OomphLibError(error_stream.str(),
1366 OOMPH_CURRENT_FUNCTION,
1367 OOMPH_EXCEPTION_LOCATION);
1368 }
1369 }
1370 }
1371
1372#endif
1373 }
1374
1375 // Reset base_mesh information
1376 Base_mesh_element_pt.clear();
1377 Base_mesh_element_pt.resize(nel_base_new, 0);
1379
1380 // Now enumerate the new base/root elements consistently
1381 unsigned count = 0;
1382 for (unsigned e = 0; e < nel_base_old; e++)
1383 {
1384 // Old root is non-halo: Just add the new roots into the
1385 // new lookup scheme consecutively
1386 if (!old_root_is_halo_or_non_existent[e])
1387 {
1388 // Loop over new root/base element
1389 unsigned n_new_root =
1390 new_base_element_associated_with_old_base_element[e].size();
1391 for (unsigned j = 0; j < n_new_root; j++)
1392 {
1393 // Store new root/base element
1394 GeneralisedElement* el_pt =
1395 new_base_element_associated_with_old_base_element[e][j];
1396 Base_mesh_element_pt[count] = el_pt;
1397 Base_mesh_element_number_plus_one[el_pt] = count + 1;
1398
1399 // Bump counter
1400 count++;
1401 }
1402 }
1403 // Old root element is halo so skip insertion (i.e. leave
1404 // entries in lookup schemes nulled) but increase counter to
1405 // ensure consistency between processors
1406 else
1407 {
1408 unsigned nskip = n_new_root[e];
1409 count += nskip;
1410 }
1411 }
1412
1413 // Re-setup the map between "root" element and number in global mesh
1414 // (used in the load_balance() routines)
1416
1417
1419 {
1420 t_end = TimingHelpers::timer();
1421 oomph_info << "Time for finishing off base mesh info "
1422 << "Problem::prune_halo_elements_and_nodes(): "
1423 << t_end - t_start << std::endl;
1424 t_start = TimingHelpers::timer();
1425 }
1426
1427
1428 // Call actions after distribute
1430
1431
1433 {
1434 t_end = TimingHelpers::timer();
1435 oomph_info << "Time for actions_after_distribute() "
1436 << "Problem::prune_halo_elements_and_nodes(): "
1437 << t_end - t_start << std::endl;
1438 t_start = TimingHelpers::timer();
1439 }
1440
1441
1442 // Re-assign the equation numbers (incl synchronisation if reqd)
1443#ifdef PARANOID
1444 unsigned n_dof = assign_eqn_numbers();
1445#else
1447#endif
1448
1449
1451 {
1452 t_end = TimingHelpers::timer();
1453 oomph_info << "Time for assign_eqn_numbers() "
1454 << "Problem::prune_halo_elements_and_nodes(): "
1455 << t_end - t_start << std::endl;
1456 t_start = TimingHelpers::timer();
1457 }
1458
1459
1460#ifdef PARANOID
1462 {
1463 if (n_dof != old_ndof)
1464 {
1465 std::ostringstream error_stream;
1466 error_stream
1467 << "Number of dofs in prune_halo_elements_and_nodes() has "
1468 "changed "
1469 << "from " << old_ndof << " to " << n_dof << "\n"
1470 << "Check that you've implemented any necessary "
1471 "actions_before/after"
1472 << "\nadapt/distribute functions, e.g. to pin redundant pressure"
1473 << " dofs etc.\n";
1474 throw OomphLibError(error_stream.str(),
1475 OOMPH_CURRENT_FUNCTION,
1476 OOMPH_EXCEPTION_LOCATION);
1477 }
1478 }
1479#endif
1480 }
1481 }
1482 }
1483
1484
1485#endif
1486
1487
1488 //===================================================================
1489 /// Build a single (global) mesh from a number
1490 /// of submeshes which are passed as a vector of pointers to the
1491 /// submeshes. The ordering is not necessarily optimal.
1492 //==============================================================
1494 {
1495#ifdef PARANOID
1496 // Has a global mesh already been built
1497 if (Mesh_pt != 0)
1498 {
1499 std::string error_message = "Problem::build_global_mesh() called,\n";
1500 error_message += " but a global mesh has already been built:\n";
1501 error_message += "Problem::Mesh_pt is not zero!\n";
1502
1503 throw OomphLibError(
1504 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1505 }
1506 // Check that there are submeshes
1507 if (Sub_mesh_pt.size() == 0)
1508 {
1509 std::string error_message = "Problem::build_global_mesh() called,\n";
1510 error_message += " but there are no submeshes:\n";
1511 error_message += "Problem::Sub_mesh_pt has no entries\n";
1512
1513 throw OomphLibError(
1514 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1515 }
1516#endif
1517
1518 // Create an empty mesh
1519 Mesh_pt = new Mesh();
1520
1521 // Call the rebuild function to construct the mesh
1523 }
1524
1525 //====================================================================
1526 /// If one of the submeshes has changed (e.g. by
1527 /// mesh adaptation) we need to update the global mesh.
1528 /// \b Note: The nodes boundary information refers to the
1529 /// boundary numbers within the submesh!
1530 /// N.B. This is essentially the same function as the Mesh constructor
1531 /// that assembles a single global mesh from submeshes
1532 //=====================================================================
1534 {
1535 // Use the function in mesh to merge the submeshes into this one
1537 }
1538
1539
1540 //================================================================
1541 /// Add a timestepper to the problem. The function will automatically
1542 /// create or resize the Time object so that it contains the appropriate
1543 /// number of levels of storage.
1544 //================================================================
1545 void Problem::add_time_stepper_pt(TimeStepper* const& time_stepper_pt)
1546 {
1547 // Add the timestepper to the vector
1549
1550 // Find the number of timesteps required by the timestepper
1551 unsigned ndt = time_stepper_pt->ndt();
1552
1553 // If time has not been allocated, create time object with the
1554 // required number of time steps
1555 if (Time_pt == 0)
1556 {
1557 Time_pt = new Time(ndt);
1558 oomph_info << "Created Time with " << ndt << " timesteps" << std::endl;
1559 }
1560 else
1561 {
1562 // If the required number of time steps is greater than currently stored
1563 // resize the time storage
1564 if (ndt > Time_pt->ndt())
1565 {
1566 Time_pt->resize(ndt);
1567 oomph_info << "Resized Time to include " << ndt << " timesteps"
1568 << std::endl;
1569 }
1570 // Otherwise report that we are OK
1571 else
1572 {
1573 oomph_info << "Time object already has storage for " << ndt
1574 << " timesteps" << std::endl;
1575 }
1576 }
1577
1578 // Pass the pointer to time to the timestepper
1580 }
1581
1582 //================================================================
1583 /// Set the explicit time stepper for the problem and also
1584 /// ensure that a time object has been created.
1585 //================================================================
1587 ExplicitTimeStepper* const& explicit_time_stepper_pt)
1588 {
1589 // Set the explicit time stepper
1591
1592 // If time has not been allocated, create time object with the
1593 // required number of time steps
1594 if (Time_pt == 0)
1595 {
1596 Time_pt = new Time(0);
1597 oomph_info << "Created Time with storage for no previous timestep"
1598 << std::endl;
1599 }
1600 else
1601 {
1602 oomph_info << "Time object already exists " << std::endl;
1603 }
1604 }
1605
1606
1607#ifdef OOMPH_HAS_MPI
1608
1609 //================================================================
1610 /// Set default first and last elements for parallel assembly
1611 /// of non-distributed problem.
1612 //================================================================
1614 {
1616 {
1617 // Minimum number of elements per processor if there are fewer elements
1618 // than processors
1619 unsigned min_el = 10;
1620
1621 // Resize and make default assignments
1622 int n_proc = this->communicator_pt()->nproc();
1623 unsigned n_elements = Mesh_pt->nelement();
1624 First_el_for_assembly.resize(n_proc, 0);
1625 Last_el_plus_one_for_assembly.resize(n_proc, 0);
1626
1627 // In the absence of any better knowledge distribute work evenly
1628 // over elements
1629 unsigned range = 0;
1630 unsigned lo_proc = 0;
1631 unsigned hi_proc = n_proc - 1;
1632 if (int(n_elements) >= n_proc)
1633 {
1634 range = unsigned(double(n_elements) / double(n_proc));
1635 }
1636 else
1637 {
1638 range = min_el;
1639 lo_proc = 0;
1640 hi_proc = unsigned(double(n_elements) / double(min_el));
1641 }
1642
1643 for (int p = lo_proc; p <= int(hi_proc); p++)
1644 {
1645 First_el_for_assembly[p] = p * range;
1646
1647 unsigned last_el_plus_one = (p + 1) * range;
1648 if (last_el_plus_one > n_elements) last_el_plus_one = n_elements;
1649 Last_el_plus_one_for_assembly[p] = last_el_plus_one;
1650 }
1651
1652 // Last one needs to incorporate any dangling elements
1653 if (int(n_elements) >= n_proc)
1654 {
1655 Last_el_plus_one_for_assembly[n_proc - 1] = n_elements;
1656 }
1657
1658 // Doc
1659 if (n_proc > 1)
1660 {
1662 {
1663 oomph_info << "Problem is not distributed. Parallel assembly of "
1664 << "Jacobian uses default partitioning: " << std::endl;
1665 for (int p = 0; p < n_proc; p++)
1666 {
1668 {
1669 oomph_info << "Proc " << p << " assembles from element "
1670 << First_el_for_assembly[p] << " to "
1671 << Last_el_plus_one_for_assembly[p] - 1 << " \n";
1672 }
1673 else
1674 {
1675 oomph_info << "Proc " << p << " assembles no elements\n";
1676 }
1677 }
1678 }
1679 }
1680 }
1681 }
1682
1683
1684 //=======================================================================
1685 /// Helper function to re-assign the first and last elements to be
1686 /// assembled by each processor during parallel assembly for
1687 /// non-distributed problem.
1688 //=======================================================================
1690 {
1691 // Wait until all processes have completed/timed their assembly
1692 MPI_Barrier(this->communicator_pt()->mpi_comm());
1693
1694 // Storage for number of processors and current processor
1695 int n_proc = this->communicator_pt()->nproc();
1696 int rank = this->communicator_pt()->my_rank();
1697
1698 // Don't bother to update if we've got fewer elements than
1699 // processors
1700 unsigned nel = Elemental_assembly_time.size();
1701 if (int(nel) < n_proc)
1702 {
1703 oomph_info << "Not re-computing distribution of elemental assembly\n"
1704 << "because there are fewer elements than processors\n";
1705 return;
1706 }
1707
1708 // Setup vectors storing the number of element timings to be sent
1709 // and the offset in the final vector
1710 Vector<int> receive_count(n_proc);
1711 Vector<int> displacement(n_proc);
1712 int offset = 0;
1713 for (int p = 0; p < n_proc; p++)
1714 {
1715 // Default distribution of labour
1716 unsigned el_lo = First_el_for_assembly[p];
1717 unsigned el_hi = Last_el_plus_one_for_assembly[p] - 1;
1718
1719 // Number of timings to be sent and offset from start in
1720 // final vector
1721 receive_count[p] = el_hi - el_lo + 1;
1722 displacement[p] = offset;
1723 offset += el_hi - el_lo + 1;
1724 }
1725
1726 // Make temporary c-style array to avoid over-writing in Gatherv below
1727 double* el_ass_time = new double[nel];
1728 for (unsigned e = 0; e < nel; e++)
1729 {
1730 el_ass_time[e] = Elemental_assembly_time[e];
1731 }
1732
1733 // Gather timings on root processor
1734 unsigned nel_local =
1736 MPI_Gatherv(&el_ass_time[First_el_for_assembly[rank]],
1737 nel_local,
1738 MPI_DOUBLE,
1740 &receive_count[0],
1741 &displacement[0],
1742 MPI_DOUBLE,
1743 0,
1744 this->communicator_pt()->mpi_comm());
1745 delete[] el_ass_time;
1746
1747 // Vector of first and last elements for each processor
1748 Vector<Vector<int>> first_and_last_element(n_proc);
1749 for (int p = 0; p < n_proc; p++)
1750 {
1751 first_and_last_element[p].resize(2);
1752 }
1753
1754 // Re-distribute work
1755 if (rank == 0)
1756 {
1758 {
1760 << std::endl
1761 << "Re-assigning distribution of element assembly over processors:"
1762 << std::endl;
1763 }
1764
1765 // Get total assembly time
1766 double total = 0.0;
1767 unsigned n_elements = Mesh_pt->nelement();
1768 for (unsigned e = 0; e < n_elements; e++)
1769 {
1770 total += Elemental_assembly_time[e];
1771 }
1772
1773 // Target load per processor
1774 double target_load = total / double(n_proc);
1775
1776 // We're on the root processor: Always start with the first element
1777 int proc = 0;
1778 first_and_last_element[0][0] = 0;
1779
1780 // Highest element we can help ourselves to if we want to leave
1781 // at least one element for all subsequent processors
1782 unsigned max_el_avail = n_elements - n_proc;
1783
1784 // Initialise total work allocated
1785 total = 0.0;
1786 for (unsigned e = 0; e < n_elements; e++)
1787 {
1788 total += Elemental_assembly_time[e];
1789
1790 // Once we have reached the target load or we've used up practically
1791 // all the elements...
1792 if ((total > target_load) || (e == max_el_avail))
1793 {
1794 // Last element for current processor
1795 first_and_last_element[proc][1] = e;
1796
1797 // Provided that we are not on the last processor
1798 if (proc < (n_proc - 1))
1799 {
1800 // Set first element for next one
1801 first_and_last_element[proc + 1][0] = e + 1;
1802
1803 // Move on to the next processor
1804 proc++;
1805 }
1806
1807 // Can have one more...
1808 max_el_avail++;
1809
1810 // Re-initialise the time
1811 total = 0.0;
1812 } // end of test for "total exceeds target"
1813 }
1814
1815
1816 // Last element for last processor
1817 first_and_last_element[n_proc - 1][1] = n_elements - 1;
1818
1819
1820 // The following block should probably be paranoidified away
1821 // but we've screwed the logic up so many times that I feel
1822 // it's safer to keep it...
1823 bool wrong = false;
1824 std::ostringstream error_stream;
1825 for (int p = 0; p < n_proc - 1; p++)
1826 {
1827 unsigned first_of_current = first_and_last_element[p][0];
1828 unsigned last_of_current = first_and_last_element[p][1];
1829 if (first_of_current > last_of_current)
1830 {
1831 wrong = true;
1832 error_stream << "Error: First/last element of proc " << p << ": "
1833 << first_of_current << " " << last_of_current
1834 << std::endl;
1835 }
1836 unsigned first_of_next = first_and_last_element[p + 1][0];
1837 if (first_of_next != (last_of_current + 1))
1838 {
1839 wrong = true;
1840 error_stream << "Error: First element of proc " << p + 1 << ": "
1841 << first_of_next << " and last element of proc " << p
1842 << ": " << last_of_current << std::endl;
1843 }
1844 }
1845 if (wrong)
1846 {
1847 throw OomphLibError(
1848 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1849 }
1850
1851
1852 // THIS TIDY UP SHOULD NO LONGER BE REQUIRED AND CAN GO AT SOME POINT
1853
1854 // //If we haven't got to the end of the processor list then
1855 // //need to shift things about slightly because the processors
1856 // //at the end will be empty.
1857 // //This can occur when you have very fast assembly times and the
1858 // //rounding errors mean that the targets are achieved before all
1859 // processors
1860 // //have been visited.
1861 // //Happens a lot when you massively oversubscribe the CPUs (which was
1862 // //only ever for testing!)
1863 // if (proc!=n_proc-1)
1864 // {
1865 // oomph_info
1866 // << "First pass did not allocate elements on every processor\n";
1867 // oomph_info <<
1868 // "Moving elements so that each processor has at least one\n";
1869
1870 // //Work out number of empty processos
1871 // unsigned n_empty_processors = n_proc - proc + 1;
1872
1873 // //Loop over the processors that do have elements
1874 // //and work out how many we need to steal elements from
1875 // unsigned n_element_on_processors=0;
1876 // do
1877 // {
1878 // //Step down the processors
1879 // --proc;
1880 // //Add the current processor to the number of empty processors
1881 // //because the elements have to be shared between processors
1882 // //including the one(s) on which they are currently stored.
1883 // ++n_empty_processors;
1884 // n_element_on_processors +=
1885 // (first_and_last_element[proc][1] -
1886 // first_and_last_element[proc][0] + 1);
1887 // }
1888 // while(n_element_on_processors < n_empty_processors);
1889
1890 // //Should now be able to put one element on each processor
1891 // //Start from the end and do so
1892 // unsigned current_element = n_elements-1;
1893 // for(int p=n_proc-1;p>proc;p--)
1894 // {
1895 // first_and_last_element[p][1] = current_element;
1896 // first_and_last_element[p][0] = --current_element;
1897 // }
1898
1899 // //Now for the last processor we touched, just adjust the final value
1900 // first_and_last_element[proc][1] = current_element;
1901 // }
1902 // //Otherwise just put the rest of the elements on the final
1903 // //processor
1904 // else
1905 // {
1906 // // Last one
1907 // first_and_last_element[n_proc-1][1]=n_elements-1;
1908 // }
1909
1910
1911 // END PRESUMED-TO-BE-UNNECESSARY BLOCK...
1912
1913
1914 // Now communicate the information
1915
1916 // Set local informationt for this (root) processor
1917 First_el_for_assembly[0] = first_and_last_element[0][0];
1918 Last_el_plus_one_for_assembly[0] = first_and_last_element[0][1] + 1;
1919
1921 {
1922 oomph_info << "Processor " << 0 << " assembles Jacobians"
1923 << " from elements " << first_and_last_element[0][0]
1924 << " to " << first_and_last_element[0][1] << " "
1925 << std::endl;
1926 }
1927
1928 // Only now can we send the information to the other processors
1929 for (int p = 1; p < n_proc; ++p)
1930 {
1931 MPI_Send(&first_and_last_element[p][0],
1932 2,
1933 MPI_INT,
1934 p,
1935 0,
1936 this->communicator_pt()->mpi_comm());
1937
1938
1940 {
1941 oomph_info << "Processor " << p << " assembles Jacobians"
1942 << " from elements " << first_and_last_element[p][0]
1943 << " to " << first_and_last_element[p][1] << " "
1944 << std::endl;
1945 }
1946 }
1947 }
1948 // Receive first and last element from root on non-master processors
1949 else
1950 {
1951 Vector<int> aux(2);
1952 MPI_Status status;
1953 MPI_Recv(&aux[0],
1954 2,
1955 MPI_INT,
1956 0,
1957 0,
1958 this->communicator_pt()->mpi_comm(),
1959 &status);
1960 First_el_for_assembly[rank] = aux[0];
1961 Last_el_plus_one_for_assembly[rank] = aux[1] + 1;
1962 }
1963
1964 // Wipe all others
1965 for (int p = 0; p < n_proc; p++)
1966 {
1967 if (p != rank)
1968 {
1969 First_el_for_assembly[p] = 0;
1971 }
1972 }
1973
1974 // The equations assembled by this processor may have changed so
1975 // we must resize the sparse assemble with arrays previous allocation
1977 }
1978
1979#endif
1980
1981 //================================================================
1982 /// Assign all equation numbers for problem: Deals with global
1983 /// data (= data that isn't attached to any elements) and then
1984 /// does the equation numbering for the elements. Bool argument
1985 /// can be set to false to ignore assigning local equation numbers
1986 /// (necessary in the parallel implementation of locate_zeta
1987 /// between multiple meshes).
1988 //================================================================
1990 const bool& assign_local_eqn_numbers)
1991 {
1992 // Check that the global mesh has been build
1993#ifdef PARANOID
1994 if (Mesh_pt == 0)
1995 {
1996 std::ostringstream error_stream;
1997 error_stream << "Global mesh does not exist, so equation numbers cannot "
1998 "be assigned.\n";
1999 // Check for sub meshes
2000 if (nsub_mesh() == 0)
2001 {
2002 error_stream << "There aren't even any sub-meshes in the Problem.\n"
2003 << "You can set the global mesh directly by using\n"
2004 << "Problem::mesh_pt() = my_mesh_pt;\n"
2005 << "OR you can use Problem::add_sub_mesh(mesh_pt); "
2006 << "to add a sub mesh.\n";
2007 }
2008 else
2009 {
2010 error_stream << "There are " << nsub_mesh() << " sub-meshes.\n";
2011 }
2012 error_stream << "You need to call Problem::build_global_mesh() to create "
2013 "a global mesh\n"
2014 << "from the sub-meshes.\n\n";
2015
2016 throw OomphLibError(
2017 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2018 }
2019#endif
2020
2021 // Number of submeshes
2022 unsigned n_sub_mesh = Sub_mesh_pt.size();
2023
2024#ifdef OOMPH_HAS_MPI
2025
2026 // Storage for number of processors
2027 int n_proc = this->communicator_pt()->nproc();
2028
2029
2030 if (n_proc > 1)
2031 {
2032 // Force re-analysis of time spent on assembly each
2033 // elemental Jacobian
2036 }
2037 else
2038 {
2040 }
2041
2042 // Re-distribution of elements over processors during assembly
2043 // must be recomputed
2045 {
2046 // Set default first and last elements for parallel assembly
2047 // of non-distributed problem.
2049 }
2050
2051#endif
2052
2053
2054 double t_start = 0.0;
2056 {
2057 t_start = TimingHelpers::timer();
2058 }
2059
2060 // Loop over all elements in the mesh and set up any additional
2061 // dependencies that they may have (e.g. storing the geometric
2062 // Data, i.e. Data that affects an element's shape in elements
2063 // with algebraic node-update functions
2064 unsigned nel = Mesh_pt->nelement();
2065 for (unsigned e = 0; e < nel; e++)
2066 {
2068 }
2069
2070#ifdef OOMPH_HAS_MPI
2071 // Complete setup of dependencies for external halo elements too
2072 unsigned n_mesh = this->nsub_mesh();
2073 for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
2074 {
2075 for (int iproc = 0; iproc < n_proc; iproc++)
2076 {
2077 unsigned n_ext_halo_el = mesh_pt(i_mesh)->nexternal_halo_element(iproc);
2078 for (unsigned e = 0; e < n_ext_halo_el; e++)
2079 {
2080 mesh_pt(i_mesh)
2081 ->external_halo_element_pt(iproc, e)
2083 }
2084 }
2085 }
2086#endif
2087
2088
2089 double t_end = 0.0;
2091 {
2092 t_end = TimingHelpers::timer();
2094 << "Time for complete setup of dependencies in assign_eqn_numbers: "
2095 << t_end - t_start << std::endl;
2096 }
2097
2098
2099 // Initialise number of dofs for reserve below
2100 unsigned n_dof = 0;
2101
2102 // Potentially loop over remainder of routine, possible re-visiting all
2103 // those parts that must be redone, following the removal of duplicate
2104 // external halo data.
2105 for (unsigned loop_count = 0; loop_count < 2; loop_count++)
2106 {
2107 //(Re)-set the dof pointer to zero length because entries are
2108 // pushed back onto it -- if it's not reset here then we get into
2109 // trouble during mesh refinement when we reassign all dofs
2110 Dof_pt.resize(0);
2111
2112 // Reserve from previous allocation if we're going around again
2113 Dof_pt.reserve(n_dof);
2114
2115 // Reset the equation number
2116 unsigned long equation_number = 0;
2117
2118 // Now set equation numbers for the global Data
2119 unsigned Nglobal_data = nglobal_data();
2120 for (unsigned i = 0; i < Nglobal_data; i++)
2121 {
2122 Global_data_pt[i]->assign_eqn_numbers(equation_number, Dof_pt);
2123 }
2124
2126 {
2127 t_start = TimingHelpers::timer();
2128 }
2129
2130 // Call assign equation numbers on the global mesh
2132
2133 // Deal with the spine meshes additional numbering
2134 // If there is only one mesh
2135 if (n_sub_mesh == 0)
2136 {
2137 if (SpineMesh* const spine_mesh_pt = dynamic_cast<SpineMesh*>(Mesh_pt))
2138 {
2139 n_dof = spine_mesh_pt->assign_global_spine_eqn_numbers(Dof_pt);
2140 }
2141 }
2142 // Otherwise loop over the sub meshes
2143 else
2144 {
2145 // Assign global equation numbers first
2146 for (unsigned i = 0; i < n_sub_mesh; i++)
2147 {
2148 if (SpineMesh* const spine_mesh_pt =
2149 dynamic_cast<SpineMesh*>(Sub_mesh_pt[i]))
2150 {
2151 n_dof = spine_mesh_pt->assign_global_spine_eqn_numbers(Dof_pt);
2152 }
2153 }
2154 }
2155
2157 {
2158 t_end = TimingHelpers::timer();
2160 << "Time for assign_global_eqn_numbers in assign_eqn_numbers: "
2161 << t_end - t_start << std::endl;
2162 t_start = TimingHelpers::timer();
2163 }
2164
2165
2166#ifdef OOMPH_HAS_MPI
2167
2168 // reset previous allocation
2170
2171 // Only synchronise if the problem has actually been
2172 // distributed.
2174 {
2175 // Synchronise the equation numbers and return the total
2176 // number of degrees of freedom in the overall problem
2177 // Do not assign local equation numbers -- we're doing this
2178 // below.
2179 n_dof = synchronise_eqn_numbers(false);
2180 }
2181 // ..else just setup the Dof_distribution_pt
2182 // NOTE: this is setup by synchronise_eqn_numbers(...)
2183 // if Problem_has_been_distributed
2184 else
2185#endif
2186 {
2188 }
2189
2191 {
2192 t_end = TimingHelpers::timer();
2193 oomph_info << "Time for Problem::synchronise_eqn_numbers in "
2194 << "Problem::assign_eqn_numbers: " << t_end - t_start
2195 << std::endl;
2196 }
2197
2198
2199#ifdef OOMPH_HAS_MPI
2200
2201
2202 // Now remove duplicate data in external halo elements
2204 {
2206 {
2207 t_start = TimingHelpers::timer();
2208 }
2209
2210 // Monitor if we've actually changed anything
2211 bool actually_removed_some_data = false;
2212
2213 // Only do it once!
2214 if (loop_count == 0)
2215 {
2216 if (n_sub_mesh == 0)
2217 {
2218 remove_duplicate_data(Mesh_pt, actually_removed_some_data);
2219 }
2220 else
2221 {
2222 for (unsigned i = 0; i < n_sub_mesh; i++)
2223 {
2224 bool tmp_actually_removed_some_data = false;
2226 tmp_actually_removed_some_data);
2227 if (tmp_actually_removed_some_data)
2228 actually_removed_some_data = true;
2229 }
2230 }
2231 }
2232
2233
2235 {
2236 t_end = TimingHelpers::timer();
2237 std::stringstream tmp;
2238 tmp << "Time for calls to Problem::remove_duplicate_data in "
2239 << "Problem::assign_eqn_numbers: " << t_end - t_start
2240 << " ; have ";
2241 if (!actually_removed_some_data)
2242 {
2243 tmp << " not ";
2244 }
2245 tmp << " removed some/any data.\n";
2246 oomph_info << tmp.str();
2247 t_start = TimingHelpers::timer();
2248 }
2249
2250 // Break out of the loop if we haven't done anything here.
2251 unsigned status = 0;
2252 if (actually_removed_some_data) status = 1;
2253
2254 // Allreduce to check if anyone has removed any data
2255 unsigned overall_status = 0;
2256 MPI_Allreduce(&status,
2257 &overall_status,
2258 1,
2259 MPI_UNSIGNED,
2260 MPI_MAX,
2261 this->communicator_pt()->mpi_comm());
2262
2263
2265 {
2266 t_end = TimingHelpers::timer();
2267 std::stringstream tmp;
2268 tmp
2269 << "Time for MPI_Allreduce after Problem::remove_duplicate_data in "
2270 << "Problem::assign_eqn_numbers: " << t_end - t_start << std::endl;
2271 oomph_info << tmp.str();
2272 t_start = TimingHelpers::timer();
2273 }
2274
2275 // Bail out if we haven't done anything here
2276 if (overall_status != 1)
2277 {
2278 break;
2279 }
2280
2281 // Big tidy up: Remove null pointers from halo/haloed node storage
2282 // for all meshes (this involves comms and therefore must be
2283 // performed outside loop over meshes so the all-to-all is only
2284 // done once)
2286
2287 // Time it...
2289 {
2290 double t_end = TimingHelpers::timer();
2291 oomph_info << "Total time for "
2292 << "Problem::remove_null_pointers_from_external_halo_node_"
2293 "storage(): "
2294 << t_end - t_start << std::endl;
2295 }
2296 }
2297 else
2298 {
2299 // Problem not distributed; no need for another loop
2300 break;
2301 }
2302
2303#else
2304
2305 // Serial run: Again no need for a second loop
2306 break;
2307
2308#endif
2309
2310 } // end of loop over fcts that need to be re-executed if
2311 // we've removed duplicate external data
2312
2313
2314 // Resize the sparse assemble with arrays previous allocation
2316
2317
2319 {
2320 t_start = TimingHelpers::timer();
2321 }
2322
2323 // Finally assign local equations
2324 if (assign_local_eqn_numbers)
2325 {
2326 if (n_sub_mesh == 0)
2327 {
2329 }
2330 else
2331 {
2332 for (unsigned i = 0; i < n_sub_mesh; i++)
2333 {
2334 Sub_mesh_pt[i]->assign_local_eqn_numbers(
2336 }
2337 }
2338 }
2339
2341 {
2342 t_end = TimingHelpers::timer();
2343 oomph_info << "Total time for all Mesh::assign_local_eqn_numbers in "
2344 << "Problem::assign_eqn_numbers: " << t_end - t_start
2345 << std::endl;
2346 }
2347
2348
2349 // and return the total number of DOFs
2350 return n_dof;
2351 }
2352 //================================================================
2353 /// Function to describe the dofs in terms of the global
2354 /// equation number, i.e. what type of value (nodal value of
2355 /// a Node; value in a Data object; value of internal Data in an
2356 /// element; etc) is the unknown with a certain global equation number.
2357 /// Output stream defaults to oomph_info.
2358 //================================================================
2359 void Problem::describe_dofs(std::ostream& out) const
2360 {
2361 // Check that the global mesh has been build
2362#ifdef PARANOID
2363 if (Mesh_pt == 0)
2364 {
2365 std::ostringstream error_stream;
2366 error_stream
2367 << "Global mesh does not exist, so equation numbers cannot be found.\n";
2368 // Check for sub meshes
2369 if (nsub_mesh() == 0)
2370 {
2371 error_stream << "There aren't even any sub-meshes in the Problem.\n"
2372 << "You can set the global mesh directly by using\n"
2373 << "Problem::mesh_pt() = my_mesh_pt;\n"
2374 << "OR you can use Problem::add_sub_mesh(mesh_pt); "
2375 << "to add a sub mesh.\n";
2376 }
2377 else
2378 {
2379 error_stream << "There are " << nsub_mesh() << " sub-meshes.\n";
2380 }
2381 error_stream << "You need to call Problem::build_global_mesh() to create "
2382 "a global mesh\n"
2383 << "from the sub-meshes.\n\n";
2384
2385 throw OomphLibError(
2386 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2387 }
2388#endif
2389
2390 out
2391 << "Although this program will describe the degrees of freedom in the \n"
2392 << "problem, it will do so using the typedef for the elements. This is \n"
2393 << "not neccesarily human readable, but there is a solution.\n"
2394 << "Pipe your program's output through c++filt, with the argument -t.\n"
2395 << "e.g. \"./two_d_multi_poisson | c++filt -t > ReadableOutput.txt\".\n "
2396 << "(Disregarding the quotes)\n\n\n";
2397
2398 out << "Classifying Global Equation Numbers" << std::endl;
2399 out << std::string(80, '-') << std::endl;
2400
2401 // Number of submeshes
2402 unsigned n_sub_mesh = Sub_mesh_pt.size();
2403
2404 // Classify Global dofs
2405 unsigned Nglobal_data = nglobal_data();
2406 for (unsigned i = 0; i < Nglobal_data; i++)
2407 {
2408 std::stringstream conversion;
2409 conversion << " in Global Data " << i << ".";
2410 std::string in(conversion.str());
2411 Global_data_pt[i]->describe_dofs(out, in);
2412 }
2413
2414 // Put string in limiting scope.
2415 {
2416 // Descend into assignment for mesh.
2417 std::string in(" in Problem's Only Mesh.");
2418 Mesh_pt->describe_dofs(out, in);
2419 }
2420
2421 // Deal with the spine meshes additional numbering:
2422 // If there is only one mesh:
2423 if (n_sub_mesh == 0)
2424 {
2425 if (SpineMesh* const spine_mesh_pt = dynamic_cast<SpineMesh*>(Mesh_pt))
2426 {
2427 std::string in(" in Problem's Only SpineMesh.");
2428 spine_mesh_pt->describe_spine_dofs(out, in);
2429 }
2430 }
2431 // Otherwise loop over the sub meshes
2432 else
2433 {
2434 // Assign global equation numbers first
2435 for (unsigned i = 0; i < n_sub_mesh; i++)
2436 {
2437 if (SpineMesh* const spine_mesh_pt =
2438 dynamic_cast<SpineMesh*>(Sub_mesh_pt[i]))
2439 {
2440 std::stringstream conversion;
2441 conversion << " in Sub-SpineMesh " << i << ".";
2442 std::string in(conversion.str());
2443 spine_mesh_pt->describe_spine_dofs(out, in);
2444 } // end if.
2445 } // end for.
2446 } // end else.
2447
2448
2449 out << std::string(80, '\\') << std::endl;
2450 out << std::string(80, '\\') << std::endl;
2451 out << std::string(80, '\\') << std::endl;
2452 out << "Classifying global eqn numbers in terms of elements." << std::endl;
2453 out << std::string(80, '-') << std::endl;
2454 out << "Eqns | Source" << std::endl;
2455 out << std::string(80, '-') << std::endl;
2456
2457 if (n_sub_mesh == 0)
2458 {
2459 std::string in(" in Problem's Only Mesh.");
2460 Mesh_pt->describe_local_dofs(out, in);
2461 }
2462 else
2463 {
2464 for (unsigned i = 0; i < n_sub_mesh; i++)
2465 {
2466 std::stringstream conversion;
2467 conversion << " in Sub-Mesh " << i << ".";
2468 std::string in(conversion.str());
2469 Sub_mesh_pt[i]->describe_local_dofs(out, in);
2470 } // End for
2471 } // End else
2472 } // End problem::describe_dofs(...)
2473
2474
2475 //================================================================
2476 /// Get the vector of dofs, i.e. a vector containing the current
2477 /// values of all unknowns.
2478 //================================================================
2480 {
2481 // Find number of dofs
2482 const unsigned long n_dof = ndof();
2483
2484 // Resize the vector
2485 dofs.build(Dof_distribution_pt, 0.0);
2486
2487 // Copy dofs into vector
2488 for (unsigned long l = 0; l < n_dof; l++)
2489 {
2490 dofs[l] = *Dof_pt[l];
2491 }
2492 }
2493
2494 /// Get history values of dofs in a double vector.
2495 void Problem::get_dofs(const unsigned& t, DoubleVector& dofs) const
2496 {
2497#ifdef PARANOID
2498 if (distributed())
2499 {
2500 throw OomphLibError("Not designed for distributed problems",
2501 OOMPH_EXCEPTION_LOCATION,
2502 OOMPH_CURRENT_FUNCTION);
2503 // might work, not sure
2504 }
2505#endif
2506
2507 // Resize the vector
2508 dofs.build(Dof_distribution_pt, 0.0);
2509
2510 // First deal with global data
2511 unsigned Nglobal_data = nglobal_data();
2512 for (unsigned i = 0; i < Nglobal_data; i++)
2513 {
2514 for (unsigned j = 0, nj = Global_data_pt[i]->nvalue(); j < nj; j++)
2515 {
2516 // For each data get the equation number and copy out the value.
2517 int eqn_number = Global_data_pt[i]->eqn_number(j);
2518 if (eqn_number >= 0)
2519 {
2520 dofs[eqn_number] = Global_data_pt[i]->value(t, j);
2521 }
2522 }
2523 }
2524
2525 // Next element internal data
2526 for (unsigned i = 0, ni = mesh_pt()->nelement(); i < ni; i++)
2527 {
2528 GeneralisedElement* ele_pt = mesh_pt()->element_pt(i);
2529 for (unsigned j = 0, nj = ele_pt->ninternal_data(); j < nj; j++)
2530 {
2531 Data* d_pt = ele_pt->internal_data_pt(j);
2532 for (unsigned k = 0, nk = d_pt->nvalue(); k < nk; k++)
2533 {
2534 int eqn_number = d_pt->eqn_number(k);
2535 if (eqn_number >= 0)
2536 {
2537 dofs[eqn_number] = d_pt->value(t, k);
2538 }
2539 }
2540 }
2541 }
2542
2543 // Now the nodes
2544 for (unsigned i = 0, ni = mesh_pt()->nnode(); i < ni; i++)
2545 {
2546 Node* node_pt = mesh_pt()->node_pt(i);
2547 for (unsigned j = 0, nj = node_pt->nvalue(); j < nj; j++)
2548 {
2549 // For each node get the equation number and copy out the value.
2550 int eqn_number = node_pt->eqn_number(j);
2551 if (eqn_number >= 0)
2552 {
2553 dofs[eqn_number] = node_pt->value(t, j);
2554 }
2555 }
2556 }
2557 }
2558
2559
2560#ifdef OOMPH_HAS_MPI
2561
2562 //=======================================================================
2563 /// Private helper function to remove repeated data
2564 /// in external haloed elements associated with specified mesh.
2565 /// Bool is true if some data was removed -- this usually requires
2566 /// re-running through certain parts of the equation numbering procedure.
2567 //======================================================================
2569 bool& actually_removed_some_data)
2570 {
2571 // // // Taken out again by MH -- clutters up output
2572 // // Doc timings if required
2573 // double t_start=0.0;
2574 // if (Global_timings::Doc_comprehensive_timings)
2575 // {
2576 // t_start=TimingHelpers::timer();
2577 // }
2578
2579 int n_proc = this->communicator_pt()->nproc();
2580 int my_rank = this->communicator_pt()->my_rank();
2581
2582 // Initialise
2583 actually_removed_some_data = false;
2584
2585 // Each individual container of external halo nodes has unique
2586 // nodes/equation numbers, but there may be some duplication between
2587 // two or more different containers; the following code checks for this
2588 // and removes the duplication by overwriting any data point with an already
2589 // existing eqn number with the original data point which had the eqn no.
2590
2591 // // Storage for existing nodes, enumerated by first non-negative
2592 // // global equation number
2593 // unsigned n_dof=ndof();
2594
2595 // Note: This used to be
2596 // Vector<Node*> global_node_pt(n_dof,0);
2597 // but this is a total killer! Memory allocation is extremely
2598 // costly and only relatively few entries are used so use
2599 // map:
2600 std::map<unsigned, Node*> global_node_pt;
2601
2602 // Only do each retained node once
2603 std::map<Node*, bool> node_done;
2604
2605 // Loop over existing "normal" elements in mesh
2606 unsigned n_element = mesh_pt->nelement();
2607 for (unsigned e = 0; e < n_element; e++)
2608 {
2609 FiniteElement* el_pt =
2610 dynamic_cast<FiniteElement*>(mesh_pt->element_pt(e));
2611 if (el_pt != 0)
2612 {
2613 // Loop over nodes
2614 unsigned n_node = el_pt->nnode();
2615 for (unsigned j = 0; j < n_node; j++)
2616 {
2617 Node* nod_pt = el_pt->node_pt(j);
2618
2619 // Have we already done the node?
2620 if (!node_done[nod_pt])
2621 {
2622 node_done[nod_pt] = true;
2623
2624 // Loop over values stored at node (if any) to find
2625 // the first non-negative eqn number
2626 unsigned first_non_negative_eqn_number_plus_one = 0;
2627 unsigned n_val = nod_pt->nvalue();
2628 for (unsigned i_val = 0; i_val < n_val; i_val++)
2629 {
2630 int eqn_no = nod_pt->eqn_number(i_val);
2631 if (eqn_no >= 0)
2632 {
2633 first_non_negative_eqn_number_plus_one = eqn_no + 1;
2634 break;
2635 }
2636 }
2637
2638 // If we haven't found a non-negative eqn number check
2639 // eqn numbers associated with solid data (if any)
2640 if (first_non_negative_eqn_number_plus_one == 0)
2641 {
2642 // Is it a solid node?
2643 SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
2644 if (solid_nod_pt != 0)
2645 {
2646 // Loop over values stored at node (if any) to find
2647 // the first non-negative eqn number
2648 unsigned n_val = solid_nod_pt->variable_position_pt()->nvalue();
2649 for (unsigned i_val = 0; i_val < n_val; i_val++)
2650 {
2651 int eqn_no =
2652 solid_nod_pt->variable_position_pt()->eqn_number(i_val);
2653 if (eqn_no >= 0)
2654 {
2655 first_non_negative_eqn_number_plus_one = eqn_no + 1;
2656 break;
2657 }
2658 }
2659 }
2660 }
2661
2662 // Associate node with first non negative global eqn number
2663 if (first_non_negative_eqn_number_plus_one > 0)
2664 {
2665 global_node_pt[first_non_negative_eqn_number_plus_one - 1] =
2666 nod_pt;
2667 }
2668
2669
2670 // Take into account master nodes too
2671 if (dynamic_cast<RefineableElement*>(el_pt) != 0)
2672 {
2673 int n_cont_int_values = dynamic_cast<RefineableElement*>(el_pt)
2674 ->ncont_interpolated_values();
2675 for (int i_cont = -1; i_cont < n_cont_int_values; i_cont++)
2676 {
2677 if (nod_pt->is_hanging(i_cont))
2678 {
2679 HangInfo* hang_pt = nod_pt->hanging_pt(i_cont);
2680 unsigned n_master = hang_pt->nmaster();
2681 for (unsigned m = 0; m < n_master; m++)
2682 {
2683 Node* master_nod_pt = hang_pt->master_node_pt(m);
2684 if (!node_done[master_nod_pt])
2685 {
2686 node_done[master_nod_pt] = true;
2687
2688 // Loop over values stored at node (if any) to find
2689 // the first non-negative eqn number
2690 unsigned first_non_negative_eqn_number_plus_one = 0;
2691 unsigned n_val = master_nod_pt->nvalue();
2692 for (unsigned i_val = 0; i_val < n_val; i_val++)
2693 {
2694 int eqn_no = master_nod_pt->eqn_number(i_val);
2695 if (eqn_no >= 0)
2696 {
2697 first_non_negative_eqn_number_plus_one = eqn_no + 1;
2698 break;
2699 }
2700 }
2701
2702 // If we haven't found a non-negative eqn number check
2703 // eqn numbers associated with solid data (if any)
2704 if (first_non_negative_eqn_number_plus_one == 0)
2705 {
2706 // If this master is a SolidNode then add its extra
2707 // eqn numbers
2708 SolidNode* master_solid_nod_pt =
2709 dynamic_cast<SolidNode*>(master_nod_pt);
2710 if (master_solid_nod_pt != 0)
2711 {
2712 // Loop over values stored at node (if any) to find
2713 // the first non-negative eqn number
2714 unsigned n_val =
2715 master_solid_nod_pt->variable_position_pt()
2716 ->nvalue();
2717 for (unsigned i_val = 0; i_val < n_val; i_val++)
2718 {
2719 int eqn_no =
2720 master_solid_nod_pt->variable_position_pt()
2721 ->eqn_number(i_val);
2722 if (eqn_no >= 0)
2723 {
2724 first_non_negative_eqn_number_plus_one =
2725 eqn_no + 1;
2726 break;
2727 }
2728 }
2729 }
2730 }
2731 // Associate node with first non negative global
2732 // eqn number
2733 if (first_non_negative_eqn_number_plus_one > 0)
2734 {
2735 global_node_pt[first_non_negative_eqn_number_plus_one -
2736 1] = master_nod_pt;
2737 }
2738
2739 } // End of not-yet-done hang node
2740 }
2741 }
2742 }
2743 }
2744 } // endif for node already done
2745 } // End of loop over nodes
2746 } // End of FiniteElement
2747
2748 // Internal data equation numbers do not need to be added since
2749 // internal data cannot be shared between distinct elements, so
2750 // internal data on locally-stored elements can never be halo.
2751 }
2752
2753 // Set to record duplicate nodes scheduled to be killed
2754 std::set<Node*> killed_nodes;
2755
2756 // Now loop over the other processors from highest to lowest
2757 // (i.e. if there is a duplicate between these containers
2758 // then this will use the node on the highest numbered processor)
2759 for (int iproc = n_proc - 1; iproc >= 0; iproc--)
2760 {
2761 // Don't have external halo elements with yourself!
2762 if (iproc != my_rank)
2763 {
2764 // Loop over external halo elements with iproc
2765 // to remove the duplicates
2766 unsigned n_element = mesh_pt->nexternal_halo_element(iproc);
2767 for (unsigned e_ext = 0; e_ext < n_element; e_ext++)
2768 {
2769 FiniteElement* finite_ext_el_pt = dynamic_cast<FiniteElement*>(
2770 mesh_pt->external_halo_element_pt(iproc, e_ext));
2771 if (finite_ext_el_pt != 0)
2772 {
2773 // Loop over nodes
2774 unsigned n_node = finite_ext_el_pt->nnode();
2775 for (unsigned j = 0; j < n_node; j++)
2776 {
2777 Node* nod_pt = finite_ext_el_pt->node_pt(j);
2778
2779 // Loop over values stored at node (if any) to find
2780 // the first non-negative eqn number
2781 unsigned first_non_negative_eqn_number_plus_one = 0;
2782 unsigned n_val = nod_pt->nvalue();
2783 for (unsigned i_val = 0; i_val < n_val; i_val++)
2784 {
2785 int eqn_no = nod_pt->eqn_number(i_val);
2786 if (eqn_no >= 0)
2787 {
2788 first_non_negative_eqn_number_plus_one = eqn_no + 1;
2789 break;
2790 }
2791 }
2792
2793 // If we haven't found a non-negative eqn number check
2794 // eqn numbers associated with solid data (if any)
2795 if (first_non_negative_eqn_number_plus_one == 0)
2796 {
2797 // Is it a solid node?
2798 SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
2799 if (solid_nod_pt != 0)
2800 {
2801 // Loop over values stored at node (if any) to find
2802 // the first non-negative eqn number
2803 unsigned n_val =
2804 solid_nod_pt->variable_position_pt()->nvalue();
2805 for (unsigned i_val = 0; i_val < n_val; i_val++)
2806 {
2807 int eqn_no =
2808 solid_nod_pt->variable_position_pt()->eqn_number(i_val);
2809 if (eqn_no >= 0)
2810 {
2811 first_non_negative_eqn_number_plus_one = eqn_no + 1;
2812 break;
2813 }
2814 }
2815 }
2816 }
2817
2818 // Identified which node we're dealing with via first non-negative
2819 // global eqn number (if there is none, everything is pinned
2820 // and we don't give a damn...)
2821 if (first_non_negative_eqn_number_plus_one > 0)
2822 {
2823 Node* existing_node_pt =
2824 global_node_pt[first_non_negative_eqn_number_plus_one - 1];
2825
2826 // Does this node already exist?
2827 if (existing_node_pt != 0)
2828 {
2829 // Record that we're about to cull one
2830 actually_removed_some_data = true;
2831
2832 // It's a duplicate, so store the duplicated one for
2833 // later killing...
2834 Node* duplicated_node_pt = nod_pt;
2835 if (!node_done[duplicated_node_pt])
2836 {
2837 // Remove node from all boundaries
2838 std::set<unsigned>* boundaries_pt;
2839 duplicated_node_pt->get_boundaries_pt(boundaries_pt);
2840 if (boundaries_pt != 0)
2841 {
2842 Vector<unsigned> bound;
2843 unsigned nb = (*boundaries_pt).size();
2844 bound.reserve(nb);
2845 for (std::set<unsigned>::iterator it =
2846 (*boundaries_pt).begin();
2847 it != (*boundaries_pt).end();
2848 it++)
2849 {
2850 bound.push_back((*it));
2851 }
2852 for (unsigned i = 0; i < nb; i++)
2853 {
2855 duplicated_node_pt);
2856 }
2857 }
2858
2859 // Get ready to kill it
2860 killed_nodes.insert(duplicated_node_pt);
2861 unsigned i_proc = unsigned(iproc);
2863 duplicated_node_pt);
2864 }
2865
2866
2867 // Note: For now we're leaving the "dangling" (no longer
2868 // accessed masters where they are; they get cleaned
2869 // up next time we delete all the external storage
2870 // for the meshes so it's a temporary "leak" only...
2871 // At some point we should probably delete them properly too
2872
2873#ifdef PARANOID
2874
2875 // Check that hang status of exiting and replacement node
2876 // matches
2877 if (dynamic_cast<RefineableElement*>(finite_ext_el_pt) != 0)
2878 {
2879 int n_cont_inter_values =
2880 dynamic_cast<RefineableElement*>(finite_ext_el_pt)
2881 ->ncont_interpolated_values();
2882 for (int i_cont = -1; i_cont < n_cont_inter_values;
2883 i_cont++)
2884 {
2885 unsigned n_master_orig = 0;
2886 if (finite_ext_el_pt->node_pt(j)->is_hanging(i_cont))
2887 {
2888 n_master_orig = finite_ext_el_pt->node_pt(j)
2889 ->hanging_pt(i_cont)
2890 ->nmaster();
2891
2892 // Temporary leak: Resolve like this:
2893 // loop over all external halo nodes and identify the
2894 // the ones that are still reached by any of the
2895 // external elements. Kill the dangling ones.
2896 }
2897 unsigned n_master_replace = 0;
2898 if (existing_node_pt->is_hanging(i_cont))
2899 {
2900 n_master_replace =
2901 existing_node_pt->hanging_pt(i_cont)->nmaster();
2902 }
2903
2904 if (n_master_orig != n_master_replace)
2905 {
2906 std::ostringstream error_stream;
2907 error_stream
2908 << "Number of master nodes for node to be replaced, "
2909 << n_master_orig << ", doesn't match"
2910 << "those of replacement node, " << n_master_replace
2911 << " for i_cont=" << i_cont << std::endl;
2912 {
2913 error_stream
2914 << "Nodal coordinates of replacement node:";
2915 unsigned ndim = existing_node_pt->ndim();
2916 for (unsigned i = 0; i < ndim; i++)
2917 {
2918 error_stream << existing_node_pt->x(i) << " ";
2919 }
2920 error_stream << "\n";
2921 error_stream << "The coordinates of its "
2922 << n_master_replace
2923 << " master nodes are: \n";
2924 for (unsigned k = 0; k < n_master_replace; k++)
2925 {
2926 Node* master_nod_pt =
2927 existing_node_pt->hanging_pt(i_cont)
2928 ->master_node_pt(k);
2929 unsigned ndim = master_nod_pt->ndim();
2930 for (unsigned i = 0; i < ndim; i++)
2931 {
2932 error_stream << master_nod_pt->x(i) << " ";
2933 }
2934 error_stream << "\n";
2935 }
2936 }
2937
2938 {
2939 error_stream
2940 << "Nodal coordinates of node to be replaced:";
2941 unsigned ndim = finite_ext_el_pt->node_pt(j)->ndim();
2942 for (unsigned i = 0; i < ndim; i++)
2943 {
2944 error_stream << finite_ext_el_pt->node_pt(j)->x(i)
2945 << " ";
2946 }
2947 error_stream << "\n";
2948 error_stream << "The coordinates of its "
2949 << n_master_orig
2950 << " master nodes are: \n";
2951 for (unsigned k = 0; k < n_master_orig; k++)
2952 {
2953 Node* master_nod_pt = finite_ext_el_pt->node_pt(j)
2954 ->hanging_pt(i_cont)
2955 ->master_node_pt(k);
2956 unsigned ndim = master_nod_pt->ndim();
2957 for (unsigned i = 0; i < ndim; i++)
2958 {
2959 error_stream << master_nod_pt->x(i) << " ";
2960 }
2961 error_stream << "\n";
2962 }
2963 }
2964
2965
2966 throw OomphLibError(error_stream.str(),
2967 OOMPH_CURRENT_FUNCTION,
2968 OOMPH_EXCEPTION_LOCATION);
2969 }
2970 }
2971 }
2972#endif
2973 // ...and point to the existing one
2974 finite_ext_el_pt->node_pt(j) = existing_node_pt;
2975 }
2976 // If it doesn't add it to the list of existing ones
2977 else
2978 {
2979 global_node_pt[first_non_negative_eqn_number_plus_one - 1] =
2980 nod_pt;
2981 node_done[nod_pt] = true;
2982 }
2983 }
2984
2985
2986 // Do the same for any master nodes of that (possibly replaced)
2987 // node
2988 if (dynamic_cast<RefineableElement*>(finite_ext_el_pt) != 0)
2989 {
2990 int n_cont_inter_values =
2991 dynamic_cast<RefineableElement*>(finite_ext_el_pt)
2992 ->ncont_interpolated_values();
2993 for (int i_cont = -1; i_cont < n_cont_inter_values; i_cont++)
2994 {
2995 if (finite_ext_el_pt->node_pt(j)->is_hanging(i_cont))
2996 {
2997 HangInfo* hang_pt =
2998 finite_ext_el_pt->node_pt(j)->hanging_pt(i_cont);
2999 unsigned n_master = hang_pt->nmaster();
3000 for (unsigned m = 0; m < n_master; m++)
3001 {
3002 Node* master_nod_pt = hang_pt->master_node_pt(m);
3003 unsigned n_val = master_nod_pt->nvalue();
3004 unsigned first_non_negative_eqn_number_plus_one = 0;
3005 for (unsigned i_val = 0; i_val < n_val; i_val++)
3006 {
3007 int eqn_no = master_nod_pt->eqn_number(i_val);
3008 if (eqn_no >= 0)
3009 {
3010 first_non_negative_eqn_number_plus_one = eqn_no + 1;
3011 break;
3012 }
3013 }
3014
3015 // If we haven't found a non-negative eqn number check
3016 // eqn numbers associated with solid data (if any)
3017 if (first_non_negative_eqn_number_plus_one == 0)
3018 {
3019 SolidNode* solid_master_nod_pt =
3020 dynamic_cast<SolidNode*>(master_nod_pt);
3021 if (solid_master_nod_pt != 0)
3022 {
3023 // Loop over values stored at node (if any) to find
3024 // the first non-negative eqn number
3025 unsigned n_val =
3026 solid_master_nod_pt->variable_position_pt()
3027 ->nvalue();
3028 for (unsigned i_val = 0; i_val < n_val; i_val++)
3029 {
3030 int eqn_no =
3031 solid_master_nod_pt->variable_position_pt()
3032 ->eqn_number(i_val);
3033 if (eqn_no >= 0)
3034 {
3035 first_non_negative_eqn_number_plus_one =
3036 eqn_no + 1;
3037 break;
3038 }
3039 }
3040 }
3041 }
3042
3043 // Identified which node we're dealing with via
3044 // first non-negative global eqn number (if there
3045 // is none, everything is pinned and we don't give a
3046 // damn...)
3047 if (first_non_negative_eqn_number_plus_one > 0)
3048 {
3049 Node* existing_node_pt = global_node_pt
3050 [first_non_negative_eqn_number_plus_one - 1];
3051
3052 // Does this node already exist?
3053 if (existing_node_pt != 0)
3054 {
3055 // Record that we're about to cull one
3056 actually_removed_some_data = true;
3057
3058 // It's a duplicate, so store the duplicated one for
3059 // later killing...
3060 Node* duplicated_node_pt = master_nod_pt;
3061
3062 if (!node_done[duplicated_node_pt])
3063 {
3064 // Remove node from all boundaries
3065 std::set<unsigned>* boundaries_pt;
3066 duplicated_node_pt->get_boundaries_pt(
3067 boundaries_pt);
3068 if (boundaries_pt != 0)
3069 {
3070 for (std::set<unsigned>::iterator it =
3071 (*boundaries_pt).begin();
3072 it != (*boundaries_pt).end();
3073 it++)
3074 {
3076 (*it), duplicated_node_pt);
3077 }
3078 }
3079
3080 killed_nodes.insert(duplicated_node_pt);
3081 unsigned i_proc = unsigned(iproc);
3083 i_proc, duplicated_node_pt);
3084 }
3085
3086 // Weight of the original node
3087 double m_weight = hang_pt->master_weight(m);
3088
3089
3090#ifdef PARANOID
3091 // Sanity check: setting replacement master
3092 // node for non-hanging node? Sign of really
3093 // f***ed up code.
3094 Node* tmp_nod_pt = finite_ext_el_pt->node_pt(j);
3095 if (!tmp_nod_pt->is_hanging(i_cont))
3096 {
3097 std::ostringstream error_stream;
3098 error_stream
3099 << "About to re-set master for i_cont= " << i_cont
3100 << " for external node (with proc " << iproc
3101 << " )" << tmp_nod_pt << " at ";
3102 unsigned n = tmp_nod_pt->ndim();
3103 for (unsigned jj = 0; jj < n; jj++)
3104 {
3105 error_stream << tmp_nod_pt->x(jj) << " ";
3106 }
3107 error_stream
3108 << " which is not hanging --> About to die!"
3109 << "Outputting offending element into oomph-info "
3110 << "stream. \n\n";
3111 oomph_info << "\n\n";
3112 finite_ext_el_pt->output(*(oomph_info.stream_pt()));
3113 oomph_info << "\n\n";
3114 oomph_info.stream_pt()->flush();
3115 throw OomphLibError(error_stream.str(),
3116 OOMPH_CURRENT_FUNCTION,
3117 OOMPH_EXCEPTION_LOCATION);
3118 }
3119#endif
3120
3121
3122 // And re-set master
3123 finite_ext_el_pt->node_pt(j)
3124 ->hanging_pt(i_cont)
3125 ->set_master_node_pt(m, existing_node_pt, m_weight);
3126 }
3127 // If it doesn't, add it to the list of existing ones
3128 else
3129 {
3130 global_node_pt
3131 [first_non_negative_eqn_number_plus_one - 1] =
3132 master_nod_pt;
3133 node_done[master_nod_pt] = true;
3134 }
3135 }
3136 } // End of loop over master nodes
3137 } // end of hanging
3138 } // end of loop over continously interpolated variables
3139 } // end refineable element (with potentially hanging node
3140
3141 } // end loop over nodes on external halo elements
3142
3143 } // End of check for finite element
3144
3145 } // end loop over external halo elements
3146 }
3147 } // end loop over processors
3148
3149
3150 // Now kill all the deleted nodes
3151 for (std::set<Node*>::iterator it = killed_nodes.begin();
3152 it != killed_nodes.end();
3153 it++)
3154 {
3155 delete (*it);
3156 }
3157
3158
3159 // oomph_info << "Number of nonzero entries in global_node_pt: "
3160 // << global_node_pt.size() << std::endl;
3161
3162 // // Time it...
3163 // // Taken out again by MH -- clutters up output
3164 // if (Global_timings::Doc_comprehensive_timings)
3165 // {
3166 // double t_end = TimingHelpers::timer();
3167 // oomph_info
3168 // << "Total time for Problem::remove_duplicate_data: "
3169 // << t_end-t_start << std::endl;
3170 // }
3171 }
3172
3173
3174 //========================================================================
3175 /// Consolidate external halo node storage by removing nulled out
3176 /// pointers in external halo and haloed schemes for all meshes.
3177 //========================================================================
3179 {
3180 // Do we have submeshes?
3181 unsigned n_mesh_loop = 1;
3182 unsigned nmesh = nsub_mesh();
3183 if (nmesh > 0)
3184 {
3185 n_mesh_loop = nmesh;
3186 }
3187
3188 // Storage for number of processors and current processor
3189 int n_proc = this->communicator_pt()->nproc();
3190 int my_rank = this->communicator_pt()->my_rank();
3191
3192 // If only one processor then return
3193 if (n_proc == 1)
3194 {
3195 return;
3196 }
3197
3198 // Loop over all (other) processors and store index of any nulled-out
3199 // external halo nodes in storage scheme.
3200
3201 // Data to be sent to each processor
3202 Vector<int> send_n(n_proc, 0);
3203
3204 // Storage for all values to be sent to all processors
3205 Vector<int> send_data;
3206
3207 // Start location within send_data for data to be sent to each processor
3208 Vector<int> send_displacement(n_proc, 0);
3209
3210 // Check missing ones
3211 for (int domain = 0; domain < n_proc; domain++)
3212 {
3213 // Set the offset for the current processor
3214 send_displacement[domain] = send_data.size();
3215
3216 // Don't bother to do anything if the processor in the loop is the
3217 // current processor
3218 if (domain != my_rank)
3219 {
3220 // Deal with sub-meshes one-by-one if required
3221 Mesh* my_mesh_pt = 0;
3222
3223 // Loop over submeshes
3224 for (unsigned imesh = 0; imesh < n_mesh_loop; imesh++)
3225 {
3226 if (nmesh == 0)
3227 {
3228 my_mesh_pt = mesh_pt();
3229 }
3230 else
3231 {
3232 my_mesh_pt = mesh_pt(imesh);
3233 }
3234
3235 // Make backup of external halo node pointers with this domain
3236 Vector<Node*> backup_pt(my_mesh_pt->external_halo_node_pt(domain));
3237
3238 // How many do we have currently?
3239 unsigned nnod = backup_pt.size();
3240
3241 // Prepare storage for updated halo nodes
3242 Vector<Node*> new_external_halo_node_pt;
3243 new_external_halo_node_pt.reserve(nnod);
3244
3245 // Loop over external halo nodes with this domain
3246 for (unsigned j = 0; j < nnod; j++)
3247 {
3248 // Get pointer to node
3249 Node* nod_pt = backup_pt[j];
3250
3251 // Has it been nulled out?
3252 if (nod_pt == 0)
3253 {
3254 // Save index of nulled out one
3255 send_data.push_back(j);
3256 }
3257 else
3258 {
3259 // Still alive: Copy across
3260 new_external_halo_node_pt.push_back(nod_pt);
3261 }
3262 }
3263
3264 // Set new external halo node vector
3265 my_mesh_pt->set_external_halo_node_pt(domain,
3266 new_external_halo_node_pt);
3267
3268 // End of data for this mesh
3269 send_data.push_back(-1);
3270
3271 } // end of loop over meshes
3272
3273 } // end skip own domain
3274
3275 // Find the number of data added to the vector
3276 send_n[domain] = send_data.size() - send_displacement[domain];
3277 }
3278
3279 // Storage for the number of data to be received from each processor
3280 Vector<int> receive_n(n_proc, 0);
3281
3282 // Now send numbers of data to be sent between all processors
3283 MPI_Alltoall(&send_n[0],
3284 1,
3285 MPI_INT,
3286 &receive_n[0],
3287 1,
3288 MPI_INT,
3289 this->communicator_pt()->mpi_comm());
3290
3291
3292 // We now prepare the data to be received
3293 // by working out the displacements from the received data
3294 Vector<int> receive_displacement(n_proc, 0);
3295 int receive_data_count = 0;
3296 for (int rank = 0; rank < n_proc; ++rank)
3297 {
3298 // Displacement is number of data received so far
3299 receive_displacement[rank] = receive_data_count;
3300 receive_data_count += receive_n[rank];
3301 }
3302
3303 // Now resize the receive buffer for all data from all processors
3304 // Make sure that it has a size of at least one
3305 if (receive_data_count == 0)
3306 {
3307 ++receive_data_count;
3308 }
3309 Vector<int> receive_data(receive_data_count);
3310
3311 // Make sure that the send buffer has size at least one
3312 // so that we don't get a segmentation fault
3313 if (send_data.size() == 0)
3314 {
3315 send_data.resize(1);
3316 }
3317
3318 // Now send the data between all the processors
3319 MPI_Alltoallv(&send_data[0],
3320 &send_n[0],
3321 &send_displacement[0],
3322 MPI_INT,
3323 &receive_data[0],
3324 &receive_n[0],
3325 &receive_displacement[0],
3326 MPI_INT,
3327 this->communicator_pt()->mpi_comm());
3328
3329 // Now use the received data
3330 for (int send_rank = 0; send_rank < n_proc; send_rank++)
3331 {
3332 // Don't bother to do anything for the processor corresponding to the
3333 // current processor or if no data were received from this processor
3334 if ((send_rank != my_rank) && (receive_n[send_rank] != 0))
3335 {
3336 // Counter for the data within the large array
3337 unsigned count = receive_displacement[send_rank];
3338
3339 // Deal with sub-meshes one-by-one if required
3340 Mesh* my_mesh_pt = 0;
3341
3342 // Loop over submeshes
3343 for (unsigned imesh = 0; imesh < n_mesh_loop; imesh++)
3344 {
3345 if (nmesh == 0)
3346 {
3347 my_mesh_pt = mesh_pt();
3348 }
3349 else
3350 {
3351 my_mesh_pt = mesh_pt(imesh);
3352 }
3353
3354 // Make backup of external haloed node pointers with this domain
3355 Vector<Node*> backup_pt =
3356 my_mesh_pt->external_haloed_node_pt(send_rank);
3357
3358 // Unpack until we reach "end of data" indicator (-1) for this mesh
3359 while (true)
3360 {
3361 // Read next entry
3362 int next_one = receive_data[count++];
3363
3364 if (next_one == -1)
3365 {
3366 break;
3367 }
3368 else
3369 {
3370 // Null out the entry
3371 backup_pt[next_one] = 0;
3372 }
3373 }
3374
3375 // How many do we have currently?
3376 unsigned nnod = backup_pt.size();
3377
3378 // Prepare storage for updated haloed nodes
3379 Vector<Node*> new_external_haloed_node_pt;
3380 new_external_haloed_node_pt.reserve(nnod);
3381
3382 // Loop over external haloed nodes with this domain
3383 for (unsigned j = 0; j < nnod; j++)
3384 {
3385 // Get pointer to node
3386 Node* nod_pt = backup_pt[j];
3387
3388 // Has it been nulled out?
3389 if (nod_pt != 0)
3390 {
3391 // Still alive: Copy across
3392 new_external_haloed_node_pt.push_back(nod_pt);
3393 }
3394 }
3395
3396 // Set new external haloed node vector
3397 my_mesh_pt->set_external_haloed_node_pt(send_rank,
3398 new_external_haloed_node_pt);
3399 }
3400 }
3401
3402 } // End of data is received
3403 }
3404
3405#endif
3406
3407
3408 //=======================================================================
3409 /// Function that sets the values of the dofs in the object
3410 //======================================================================
3412 {
3413 const unsigned long n_dof = this->ndof();
3414#ifdef PARANOID
3415 if (n_dof != dofs.nrow())
3416 {
3417 std::ostringstream error_stream;
3418 error_stream << "Number of degrees of freedom in vector argument "
3419 << dofs.nrow() << "\n"
3420 << "does not equal number of degrees of freedom in problem "
3421 << n_dof;
3422 throw OomphLibError(
3423 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3424 }
3425#endif
3426 for (unsigned long l = 0; l < n_dof; l++)
3427 {
3428 *Dof_pt[l] = dofs[l];
3429 }
3430 }
3431
3432 /// Set history values of dofs
3433 void Problem::set_dofs(const unsigned& t, DoubleVector& dofs)
3434 {
3435#ifdef PARANOID
3436 if (distributed())
3437 {
3438 throw OomphLibError("Not designed for distributed problems",
3439 OOMPH_EXCEPTION_LOCATION,
3440 OOMPH_CURRENT_FUNCTION);
3441 // might work if the dofs vector is distributed in the right way...
3442 }
3443#endif
3444
3445 // First deal with global data
3446 unsigned Nglobal_data = nglobal_data();
3447 for (unsigned i = 0; i < Nglobal_data; i++)
3448 {
3449 for (unsigned j = 0, nj = Global_data_pt[i]->nvalue(); j < nj; j++)
3450 {
3451 // For each data get the equation number and copy out the value.
3452 int eqn_number = Global_data_pt[i]->eqn_number(j);
3453 if (eqn_number >= 0)
3454 {
3455 Global_data_pt[i]->set_value(t, j, dofs[eqn_number]);
3456 }
3457 }
3458 }
3459
3460 // Next element internal data
3461 for (unsigned i = 0, ni = mesh_pt()->nelement(); i < ni; i++)
3462 {
3463 GeneralisedElement* ele_pt = mesh_pt()->element_pt(i);
3464 for (unsigned j = 0, nj = ele_pt->ninternal_data(); j < nj; j++)
3465 {
3466 Data* d_pt = ele_pt->internal_data_pt(j);
3467 for (unsigned k = 0, nk = d_pt->nvalue(); k < nk; k++)
3468 {
3469 int eqn_number = d_pt->eqn_number(k);
3470 if (eqn_number >= 0)
3471 {
3472 d_pt->set_value(t, k, dofs[eqn_number]);
3473 }
3474 }
3475 }
3476 }
3477
3478 // Now the nodes
3479 for (unsigned i = 0, ni = mesh_pt()->nnode(); i < ni; i++)
3480 {
3481 Node* node_pt = mesh_pt()->node_pt(i);
3482 for (unsigned j = 0, nj = node_pt->nvalue(); j < nj; j++)
3483 {
3484 // For each node get the equation number and copy out the value.
3485 int eqn_number = node_pt->eqn_number(j);
3486 if (eqn_number >= 0)
3487 {
3488 node_pt->set_value(t, j, dofs[eqn_number]);
3489 }
3490 }
3491 }
3492 }
3493
3494
3495 /// Set history values of dofs from the type of vector stored in
3496 /// problem::Dof_pt.
3497 void Problem::set_dofs(const unsigned& t, Vector<double*>& dof_pt)
3498 {
3499#ifdef PARANOID
3500 if (distributed())
3501 {
3502 throw OomphLibError("Not implemented for distributed problems!",
3503 OOMPH_EXCEPTION_LOCATION,
3504 OOMPH_CURRENT_FUNCTION);
3505 }
3506#endif
3507
3508 // If we have any spine meshes I think there might be more degrees
3509 // of freedom there. I don't use them though so I'll let someone who
3510 // knows what they are doing handle it. --David Shepherd
3511
3512 // First deal with global data
3513 unsigned Nglobal_data = nglobal_data();
3514 for (unsigned i = 0; i < Nglobal_data; i++)
3515 {
3516 for (unsigned j = 0, nj = Global_data_pt[i]->nvalue(); j < nj; j++)
3517 {
3518 // For each data get the equation number and copy in the value.
3519 int eqn_number = Global_data_pt[i]->eqn_number(j);
3520 if (eqn_number >= 0)
3521 {
3522 Global_data_pt[i]->set_value(t, j, *(dof_pt[eqn_number]));
3523 }
3524 }
3525 }
3526
3527 // Now the mesh data
3528 // nodes
3529 for (unsigned i = 0, ni = mesh_pt()->nnode(); i < ni; i++)
3530 {
3531 Node* node_pt = mesh_pt()->node_pt(i);
3532 for (unsigned j = 0, nj = node_pt->nvalue(); j < nj; j++)
3533 {
3534 // For each node get the equation number and copy in the value.
3535 int eqn_number = node_pt->eqn_number(j);
3536 if (eqn_number >= 0)
3537 {
3538 node_pt->set_value(t, j, *(dof_pt[eqn_number]));
3539 }
3540 }
3541 }
3542
3543 // and non-nodal data inside elements
3544 for (unsigned i = 0, ni = mesh_pt()->nelement(); i < ni; i++)
3545 {
3546 GeneralisedElement* ele_pt = mesh_pt()->element_pt(i);
3547 for (unsigned j = 0, nj = ele_pt->ninternal_data(); j < nj; j++)
3548 {
3549 Data* data_pt = ele_pt->internal_data_pt(j);
3550 // For each node get the equation number and copy in the value.
3551 int eqn_number = data_pt->eqn_number(j);
3552 if (eqn_number >= 0)
3553 {
3554 data_pt->set_value(t, j, *(dof_pt[eqn_number]));
3555 }
3556 }
3557 }
3558 }
3559
3560
3561 //===================================================================
3562 /// Function that adds the values to the dofs
3563 //==================================================================
3564 void Problem::add_to_dofs(const double& lambda,
3565 const DoubleVector& increment_dofs)
3566 {
3567 const unsigned long n_dof = this->ndof();
3568 for (unsigned long l = 0; l < n_dof; l++)
3569 {
3570 *Dof_pt[l] += lambda * increment_dofs[l];
3571 }
3572 }
3573
3574
3575 //=========================================================================
3576 /// Return the residual vector multiplied by the inverse mass matrix
3577 /// Virtual so that it can be overloaded for mpi problems
3578 //=========================================================================
3580 {
3581 // This function does not make sense for assembly handlers other than the
3582 // default, so complain if we try to call it with another handler
3583
3584#ifdef PARANOID
3585 // If we are not the default, then complain
3587 {
3588 std::ostringstream error_stream;
3589 error_stream << "The function get_inverse_mass_matrix_times_residuals() "
3590 "can only be\n"
3591 << "used with the default assembly handler\n\n";
3592 throw OomphLibError(
3593 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3594 }
3595#endif
3596
3597 // Find the number of degrees of freedom in the problem
3598 const unsigned n_dof = this->ndof();
3599
3600 // Resize the vector
3601 LinearAlgebraDistribution dist(this->communicator_pt(), n_dof, false);
3602 Mres.build(&dist, 0.0);
3603
3604 // If we have discontinuous formulation
3605 // We can invert the mass matrix element by element
3607 {
3608 // Loop over the elements and get their residuals
3609 const unsigned n_element = Problem::mesh_pt()->nelement();
3610 Vector<double> element_Mres;
3611 for (unsigned e = 0; e < n_element; e++)
3612 {
3613 // Cache the element
3614 DGElement* const elem_pt =
3615 dynamic_cast<DGElement*>(Problem::mesh_pt()->element_pt(e));
3616
3617 // Find the elemental inverse mass matrix times residuals
3618 const unsigned n_el_dofs = elem_pt->ndof();
3619 elem_pt->get_inverse_mass_matrix_times_residuals(element_Mres);
3620
3621 // Add contribution to global matrix
3622 for (unsigned i = 0; i < n_el_dofs; i++)
3623 {
3624 Mres[elem_pt->eqn_number(i)] = element_Mres[i];
3625 }
3626 }
3627 }
3628 // Otherwise it's continous and we must invert the full
3629 // mass matrix via a global linear solve.
3630 else
3631 {
3632 // Now do the linear solve -- recycling Mass matrix if requested
3633 // If we already have the factorised mass matrix, then resolve
3635 {
3637 {
3638 oomph_info << "Not recomputing Mass Matrix " << std::endl;
3639 }
3640
3641 // Get the residuals
3642 DoubleVector residuals(&dist, 0.0);
3643 this->get_residuals(residuals);
3644
3645 // Resolve the linear system
3647 residuals, Mres);
3648 }
3649 // Otherwise solve for the first time
3650 else
3651 {
3652 // If we wish to reuse the mass matrix, then enable resolve
3654 {
3656 {
3657 oomph_info << "Enabling resolve in explicit timestep" << std::endl;
3658 }
3660 ->enable_resolve();
3661 }
3662
3663 // Use a custom assembly handler to assemble and invert the mass matrix
3664
3665 // Store the old assembly handler
3666 AssemblyHandler* old_assembly_handler_pt = this->assembly_handler_pt();
3667 // Set the assembly handler to the explicit timestep handler
3669
3670 // Solve the linear system
3672 Mres);
3673 // The mass matrix has now been computed
3675
3676 // Delete the Explicit Timestep handler
3677 delete this->assembly_handler_pt();
3678 // Reset the assembly handler to the original handler
3679 this->assembly_handler_pt() = old_assembly_handler_pt;
3680 }
3681 }
3682 }
3683
3685 {
3686 // Loop over timesteppers: make them (temporarily) steady and store their
3687 // is_steady status.
3688 unsigned n_time_steppers = this->ntime_stepper();
3689 std::vector<bool> was_steady(n_time_steppers);
3690 for (unsigned i = 0; i < n_time_steppers; i++)
3691 {
3692 was_steady[i] = time_stepper_pt(i)->is_steady();
3694 }
3695
3696 // Calculate f using the residual/jacobian machinary.
3698
3699 // Reset the is_steady status of all timesteppers that weren't already
3700 // steady when we came in here and reset their weights
3701 for (unsigned i = 0; i < n_time_steppers; i++)
3702 {
3703 if (!was_steady[i])
3704 {
3706 }
3707 }
3708 }
3709
3710
3711 //================================================================
3712 /// Get the total residuals Vector for the problem
3713 //================================================================
3715 {
3716 // Three different cases; if MPI_Helpers::MPI_has_been_initialised=true
3717 // this means MPI_Helpers::init() has been called. This could happen on a
3718 // code compiled with MPI but run serially; in this instance the
3719 // get_residuals function still works on one processor.
3720 //
3721 // Secondly, if a code has been compiled with MPI, but MPI_Helpers::init()
3722 // has not been called, then MPI_Helpers::MPI_has_been_initialised=false
3723 // and the code calls...
3724 //
3725 // Thirdly, the serial version (compiled by all, but only run when compiled
3726 // with MPI if MPI_Helpers::MPI_has_been_initialised=false
3727
3728 // Check that the residuals has the correct number of rows if it has been
3729 // setup
3730#ifdef PARANOID
3731 if (residuals.built())
3732 {
3733 if (residuals.distribution_pt()->nrow() != this->ndof())
3734 {
3735 std::ostringstream error_stream;
3736 error_stream << "The distribution of the residuals vector does not "
3737 "have the correct\n"
3738 << "number of global rows\n";
3739
3740 throw OomphLibError(
3741 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3742 }
3743 }
3744#endif
3745
3746 // Find the number of rows
3747 const unsigned nrow = this->ndof();
3748
3749 // Determine the distribution for the residuals vector
3750 // IF the vector has distribution setup then use that
3751 // ELSE determine the distribution based on the
3752 // distributed_matrix_distribution enum
3753 LinearAlgebraDistribution* dist_pt = 0;
3754 if (residuals.built())
3755 {
3756 dist_pt = new LinearAlgebraDistribution(residuals.distribution_pt());
3757 }
3758 else
3759 {
3760#ifdef OOMPH_HAS_MPI
3761 // if problem is only one one processor assemble non-distributed
3762 // distribution
3763 if (Communicator_pt->nproc() == 1)
3764 {
3765 dist_pt = new LinearAlgebraDistribution(Communicator_pt, nrow, false);
3766 }
3767 // if the problem is not distributed then assemble the jacobian with
3768 // a uniform distributed distribution
3770 {
3771 dist_pt = new LinearAlgebraDistribution(Communicator_pt, nrow, true);
3772 }
3773 // otherwise the problem is a distributed problem
3774 else
3775 {
3777 {
3779 dist_pt =
3781 break;
3784 break;
3786 LinearAlgebraDistribution* uniform_dist_pt =
3788 bool use_problem_dist = true;
3789 unsigned nproc = Communicator_pt->nproc();
3790 for (unsigned p = 0; p < nproc; p++)
3791 {
3792 if ((double)Dof_distribution_pt->nrow_local(p) >
3793 ((double)uniform_dist_pt->nrow_local(p)) * 1.1)
3794 {
3795 use_problem_dist = false;
3796 }
3797 }
3798 if (use_problem_dist)
3799 {
3801 }
3802 else
3803 {
3804 dist_pt = new LinearAlgebraDistribution(uniform_dist_pt);
3805 }
3806 delete uniform_dist_pt;
3807 break;
3808 }
3809 }
3810#else
3811 dist_pt = new LinearAlgebraDistribution(Communicator_pt, nrow, false);
3812#endif
3813 }
3814
3815 // Locally cache pointer to assembly handler
3817
3818 // Build and zero the residuals
3819 residuals.build(dist_pt, 0.0);
3820
3821 // Serial (or one processor case)
3822#ifdef OOMPH_HAS_MPI
3823 if (this->communicator_pt()->nproc() == 1)
3824 {
3825#endif // OOMPH_HAS_MPI
3826 // Loop over all the elements
3827 unsigned long Element_pt_range = Mesh_pt->nelement();
3828 for (unsigned long e = 0; e < Element_pt_range; e++)
3829 {
3830 // Get the pointer to the element
3832 // Find number of dofs in the element
3833 unsigned n_element_dofs = assembly_handler_pt->ndof(elem_pt);
3834 // Set up an array
3835 Vector<double> element_residuals(n_element_dofs);
3836 // Fill the array
3837 assembly_handler_pt->get_residuals(elem_pt, element_residuals);
3838 // Now loop over the dofs and assign values to global Vector
3839 for (unsigned l = 0; l < n_element_dofs; l++)
3840 {
3841 residuals[assembly_handler_pt->eqn_number(elem_pt, l)] +=
3842 element_residuals[l];
3843 }
3844 }
3845 // Otherwise parallel case
3846#ifdef OOMPH_HAS_MPI
3847 }
3848 else
3849 {
3850 // Store the current assembly handler
3851 AssemblyHandler* const old_assembly_handler_pt = Assembly_handler_pt;
3852 // Create a new assembly handler that only assembles the residuals
3854 new ParallelResidualsHandler(old_assembly_handler_pt);
3855
3856 // Setup memory for parallel sparse assemble
3857 // No matrix so all size zero
3858 Vector<int*> column_index;
3859 Vector<int*> row_start;
3860 Vector<double*> value;
3861 Vector<unsigned> nnz;
3862 // One set of residuals of sizer one
3863 Vector<double*> res(1);
3864
3865 // Call the parallel sparse assemble, that should only assemble residuals
3867 dist_pt, column_index, row_start, value, nnz, res);
3868 // Fill in the residuals data
3869 residuals.set_external_values(res[0], true);
3870
3871 // Delete new assembly handler
3872 delete Assembly_handler_pt;
3873 // Reset the assembly handler to the original
3874 Assembly_handler_pt = old_assembly_handler_pt;
3875 }
3876#endif
3877
3878 // Delete the distribution
3879 delete dist_pt;
3880 }
3881
3882 //=============================================================================
3883 /// Get the fully assembled residual vector and Jacobian matrix
3884 /// in dense storage. The DoubleVector residuals returned will be
3885 /// non-distributed. If on calling this method the DoubleVector residuals is
3886 /// setup then it must be non-distributed and of the correct length.
3887 /// The matrix type DenseDoubleMatrix is not distributable and therefore
3888 /// the residual vector is also assumed to be non distributable.
3889 //=============================================================================
3891 DenseDoubleMatrix& jacobian)
3892 {
3893 // get the number of degrees of freedom
3894 unsigned n_dof = ndof();
3895
3896#ifdef PARANOID
3897 // PARANOID checks : if the distribution of residuals is setup then it must
3898 // must not be distributed, have the right number of rows, and the same
3899 // communicator as the problem
3900 if (residuals.built())
3901 {
3902 if (residuals.distribution_pt()->distributed())
3903 {
3904 std::ostringstream error_stream;
3905 error_stream
3906 << "If the DoubleVector residuals is setup then it must not "
3907 << "be distributed.";
3908 throw OomphLibError(
3909 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3910 }
3911 if (residuals.distribution_pt()->nrow() != n_dof)
3912 {
3913 std::ostringstream error_stream;
3914 error_stream
3915 << "If the DoubleVector residuals is setup then it must have"
3916 << " the correct number of rows";
3917 throw OomphLibError(
3918 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3919 }
3920 if (!(*Communicator_pt ==
3921 *residuals.distribution_pt()->communicator_pt()))
3922 {
3923 std::ostringstream error_stream;
3924 error_stream
3925 << "If the DoubleVector residuals is setup then it must have"
3926 << " the same communicator as the problem.";
3927 throw OomphLibError(
3928 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3929 }
3930 }
3931#endif
3932
3933 // set the residuals distribution if it is not setup
3934 if (!residuals.built())
3935 {
3936 LinearAlgebraDistribution dist(Communicator_pt, n_dof, false);
3937 residuals.build(&dist, 0.0);
3938 }
3939 // else just zero the residuals
3940 else
3941 {
3942 residuals.initialise(0.0);
3943 }
3944
3945 // Resize the matrices -- this cannot always be done externally
3946 // because get_jacobian exists in many different versions for
3947 // different storage formats -- resizing a CC or CR matrix doesn't
3948 // make sense.
3949
3950 // resize the jacobian
3951 jacobian.resize(n_dof, n_dof);
3952 jacobian.initialise(0.0);
3953
3954 // Locally cache pointer to assembly handler
3956
3957 // Loop over all the elements
3958 unsigned long n_element = Mesh_pt->nelement();
3959 for (unsigned long e = 0; e < n_element; e++)
3960 {
3961 // Get the pointer to the element
3963 // Find number of dofs in the element
3964 unsigned n_element_dofs = assembly_handler_pt->ndof(elem_pt);
3965 // Set up an array
3966 Vector<double> element_residuals(n_element_dofs);
3967 // Set up a matrix
3968 DenseMatrix<double> element_jacobian(n_element_dofs);
3969 // Fill the array
3971 elem_pt, element_residuals, element_jacobian);
3972 // Now loop over the dofs and assign values to global Vector
3973 for (unsigned l = 0; l < n_element_dofs; l++)
3974 {
3975 unsigned long eqn_number = assembly_handler_pt->eqn_number(elem_pt, l);
3976 residuals[eqn_number] += element_residuals[l];
3977 for (unsigned l2 = 0; l2 < n_element_dofs; l2++)
3978 {
3979 jacobian(eqn_number, assembly_handler_pt->eqn_number(elem_pt, l2)) +=
3980 element_jacobian(l, l2);
3981 }
3982 }
3983 }
3984 }
3985
3986 //=============================================================================
3987 /// Return the fully-assembled Jacobian and residuals for the problem,
3988 /// in the case where the Jacobian matrix is in a distributable
3989 /// row compressed storage format.
3990 /// 1. If the distribution of the jacobian and residuals is setup then, they
3991 /// will be returned with that distribution.
3992 /// Note. the jacobian and residuals must have the same distribution.
3993 /// 2. If the distribution of the jacobian and residuals are not setup then
3994 /// their distribution will computed based on:
3995 /// Distributed_problem_matrix_distribution.
3996 //=============================================================================
3998 {
3999 // Three different cases; if MPI_Helpers::MPI_has_been_initialised=true
4000 // this means MPI_Helpers::setup() has been called. This could happen on a
4001 // code compiled with MPI but run serially; in this instance the
4002 // get_residuals function still works on one processor.
4003 //
4004 // Secondly, if a code has been compiled with MPI, but MPI_Helpers::setup()
4005 // has not been called, then MPI_Helpers::MPI_has_been_initialised=false
4006 // and the code calls...
4007 //
4008 // Thirdly, the serial version (compiled by all, but only run when compiled
4009 // with MPI if MPI_Helpers::MPI_has_been_initialised=false
4010 //
4011 // The only case where an MPI code cannot run serially at present
4012 // is one where the distribute function is used (i.e. METIS is called)
4013
4014 // Allocate storage for the matrix entries
4015 // The generalised Vector<Vector<>> structure is required
4016 // for the most general interface to sparse_assemble() which allows
4017 // the assembly of multiple matrices at once.
4018 Vector<int*> column_index(1);
4019 Vector<int*> row_start(1);
4020 Vector<double*> value(1);
4021 Vector<unsigned> nnz(1);
4022
4023#ifdef PARANOID
4024 // PARANOID checks that the distribution of the jacobian matches that of the
4025 // residuals (if they are setup) and that they have the right number of rows
4026 if (residuals.built() && jacobian.distribution_built())
4027 {
4028 if (!(*residuals.distribution_pt() == *jacobian.distribution_pt()))
4029 {
4030 std::ostringstream error_stream;
4031 error_stream << "The distribution of the residuals must "
4032 << "be the same as the distribution of the jacobian.";
4033 throw OomphLibError(
4034 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4035 }
4036 if (jacobian.distribution_pt()->nrow() != this->ndof())
4037 {
4038 std::ostringstream error_stream;
4039 error_stream
4040 << "The distribution of the jacobian and residuals does not"
4041 << "have the correct number of global rows.";
4042 throw OomphLibError(
4043 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4044 }
4045 }
4046 else if (residuals.built() != jacobian.distribution_built())
4047 {
4048 std::ostringstream error_stream;
4049 error_stream << "The distribution of the jacobian and residuals must "
4050 << "both be setup or both not setup";
4051 throw OomphLibError(
4052 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4053 }
4054#endif
4055
4056
4057 // Allocate generalised storage format for passing to sparse_assemble()
4058 Vector<double*> res(1);
4059
4060 // number of rows
4061 unsigned nrow = this->ndof();
4062
4063 // determine the distribution for the jacobian.
4064 // IF the jacobian has distribution setup then use that
4065 // ELSE determine the distribution based on the
4066 // distributed_matrix_distribution enum
4067 LinearAlgebraDistribution* dist_pt = 0;
4068 if (jacobian.distribution_built())
4069 {
4070 dist_pt = new LinearAlgebraDistribution(jacobian.distribution_pt());
4071 }
4072 else
4073 {
4074#ifdef OOMPH_HAS_MPI
4075 // if problem is only one one processor
4076 if (Communicator_pt->nproc() == 1)
4077 {
4078 dist_pt = new LinearAlgebraDistribution(Communicator_pt, nrow, false);
4079 }
4080 // if the problem is not distributed then assemble the jacobian with
4081 // a uniform distributed distribution
4083 {
4084 dist_pt = new LinearAlgebraDistribution(Communicator_pt, nrow, true);
4085 }
4086 // otherwise the problem is a distributed problem
4087 else
4088 {
4090 {
4092 dist_pt =
4094 break;
4097 break;
4099 LinearAlgebraDistribution* uniform_dist_pt =
4101 bool use_problem_dist = true;
4102 unsigned nproc = Communicator_pt->nproc();
4103 for (unsigned p = 0; p < nproc; p++)
4104 {
4105 if ((double)Dof_distribution_pt->nrow_local(p) >
4106 ((double)uniform_dist_pt->nrow_local(p)) * 1.1)
4107 {
4108 use_problem_dist = false;
4109 }
4110 }
4111 if (use_problem_dist)
4112 {
4114 }
4115 else
4116 {
4117 dist_pt = new LinearAlgebraDistribution(uniform_dist_pt);
4118 }
4119 delete uniform_dist_pt;
4120 break;
4121 }
4122 }
4123#else
4124 dist_pt = new LinearAlgebraDistribution(Communicator_pt, nrow, false);
4125#endif
4126 }
4127
4128
4129 // The matrix is in compressed row format
4130 bool compressed_row_flag = true;
4131
4132#ifdef OOMPH_HAS_MPI
4133 //
4134 if (Communicator_pt->nproc() == 1)
4135 {
4136#endif
4138 column_index, row_start, value, nnz, res, compressed_row_flag);
4139 jacobian.build(dist_pt);
4140 jacobian.build_without_copy(
4141 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4142 residuals.build(dist_pt, 0.0);
4143 residuals.set_external_values(res[0], true);
4144#ifdef OOMPH_HAS_MPI
4145 }
4146 else
4147 {
4148 if (dist_pt->distributed())
4149 {
4151 dist_pt, column_index, row_start, value, nnz, res);
4152 jacobian.build(dist_pt);
4153 jacobian.build_without_copy(
4154 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4155 residuals.build(dist_pt, 0.0);
4156 residuals.set_external_values(res[0], true);
4157 }
4158 else
4159 {
4160 LinearAlgebraDistribution* temp_dist_pt =
4161 new LinearAlgebraDistribution(Communicator_pt, dist_pt->nrow(), true);
4163 temp_dist_pt, column_index, row_start, value, nnz, res);
4164 jacobian.build(temp_dist_pt);
4165 jacobian.build_without_copy(
4166 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4167 jacobian.redistribute(dist_pt);
4168 residuals.build(temp_dist_pt, 0.0);
4169 residuals.set_external_values(res[0], true);
4170 residuals.redistribute(dist_pt);
4171 delete temp_dist_pt;
4172 }
4173 }
4174#endif
4175
4176 // clean up dist_pt and residuals_vector pt
4177 delete dist_pt;
4178 }
4179
4180 //=============================================================================
4181 /// Return the fully-assembled Jacobian and residuals for the problem,
4182 /// in the case when the jacobian matrix is in column-compressed storage
4183 /// format.
4184 //=============================================================================
4186 {
4187 // Three different cases; if MPI_Helpers::MPI_has_been_initialised=true
4188 // this means MPI_Helpers::setup() has been called. This could happen on a
4189 // code compiled with MPI but run serially; in this instance the
4190 // get_residuals function still works on one processor.
4191 //
4192 // Secondly, if a code has been compiled with MPI, but MPI_Helpers::setup()
4193 // has not been called, then MPI_Helpers::MPI_has_been_initialised=false
4194 // and the code calls...
4195 //
4196 // Thirdly, the serial version (compiled by all, but only run when compiled
4197 // with MPI if MPI_Helpers::MPI_has_been_5Binitialised=false
4198 //
4199 // The only case where an MPI code cannot run serially at present
4200 // is one where the distribute function is used (i.e. METIS is called)
4201
4202 // get the number of degrees of freedom
4203 unsigned n_dof = ndof();
4204
4205#ifdef PARANOID
4206 // PARANOID checks : if the distribution of residuals is setup then it must
4207 // must not be distributed, have the right number of rows, and the same
4208 // communicator as the problem
4209 if (residuals.built())
4210 {
4211 if (residuals.distribution_pt()->distributed())
4212 {
4213 std::ostringstream error_stream;
4214 error_stream
4215 << "If the DoubleVector residuals is setup then it must not "
4216 << "be distributed.";
4217 throw OomphLibError(
4218 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4219 }
4220 if (residuals.distribution_pt()->nrow() != n_dof)
4221 {
4222 std::ostringstream error_stream;
4223 error_stream
4224 << "If the DoubleVector residuals is setup then it must have"
4225 << " the correct number of rows";
4226 throw OomphLibError(
4227 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4228 }
4229 if (!(*Communicator_pt ==
4230 *residuals.distribution_pt()->communicator_pt()))
4231 {
4232 std::ostringstream error_stream;
4233 error_stream
4234 << "If the DoubleVector residuals is setup then it must have"
4235 << " the same communicator as the problem.";
4236 throw OomphLibError(
4237 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4238 }
4239 }
4240#endif
4241
4242 // Allocate storage for the matrix entries
4243 // The generalised Vector<Vector<>> structure is required
4244 // for the most general interface to sparse_assemble() which allows
4245 // the assembly of multiple matrices at once.
4246 Vector<int*> row_index(1);
4247 Vector<int*> column_start(1);
4248 Vector<double*> value(1);
4249
4250 // Allocate generalised storage format for passing to sparse_assemble()
4251 Vector<double*> res(1);
4252
4253 // allocate storage for the number of non-zeros in each matrix
4254 Vector<unsigned> nnz(1);
4255
4256 // The matrix is in compressed column format
4257 bool compressed_row_flag = false;
4258
4259 // get the distribution for the residuals
4261 if (!residuals.built())
4262 {
4263 dist_pt =
4264 new LinearAlgebraDistribution(Communicator_pt, this->ndof(), false);
4265 }
4266 else
4267 {
4268 dist_pt = new LinearAlgebraDistribution(residuals.distribution_pt());
4269 }
4270
4271#ifdef OOMPH_HAS_MPI
4272 if (communicator_pt()->nproc() == 1)
4273 {
4274#endif
4276 row_index, column_start, value, nnz, res, compressed_row_flag);
4277 jacobian.build_without_copy(
4278 value[0], row_index[0], column_start[0], nnz[0], n_dof, n_dof);
4279 residuals.build(dist_pt, 0.0);
4280 residuals.set_external_values(res[0], true);
4281#ifdef OOMPH_HAS_MPI
4282 }
4283 else
4284 {
4285 std::ostringstream error_stream;
4286 error_stream << "Cannot assemble a CCDoubleMatrix Jacobian on more "
4287 << "than one processor.";
4288 throw OomphLibError(
4289 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4290 }
4291#endif
4292
4293 // clean up
4294 delete dist_pt;
4295 }
4296
4297
4298 //===================================================================
4299 /// Set all pinned values to zero.
4300 /// Used to set boundary conditions to be homogeneous in the copy
4301 /// of the problem used in adaptive bifurcation tracking
4302 /// (ALH: TEMPORARY HACK, WILL BE FIXED)
4303 //==================================================================
4305 {
4306 // NOTE THIS DOES NOT ZERO ANY SPINE DATA, but otherwise everything else
4307 // should be zeroed
4308
4309 // Zero any pinned global Data
4310 const unsigned n_global_data = nglobal_data();
4311 for (unsigned i = 0; i < n_global_data; i++)
4312 {
4313 Data* const local_data_pt = Global_data_pt[i];
4314 const unsigned n_value = local_data_pt->nvalue();
4315 for (unsigned j = 0; j < n_value; j++)
4316 {
4317 // If the data value is pinned set the value to zero
4318 if (local_data_pt->is_pinned(j))
4319 {
4320 local_data_pt->set_value(j, 0.0);
4321 }
4322 }
4323 }
4324
4325 // Loop over the submeshes:
4326 const unsigned n_sub_mesh = Sub_mesh_pt.size();
4327 if (n_sub_mesh == 0)
4328 {
4329 // Loop over the nodes in the element
4330 const unsigned n_node = Mesh_pt->nnode();
4331 for (unsigned n = 0; n < n_node; n++)
4332 {
4333 Node* const local_node_pt = Mesh_pt->node_pt(n);
4334 const unsigned n_value = local_node_pt->nvalue();
4335 for (unsigned j = 0; j < n_value; j++)
4336 {
4337 // If the data value is pinned set the value to zero
4338 if (local_node_pt->is_pinned(j))
4339 {
4340 local_node_pt->set_value(j, 0.0);
4341 }
4342 }
4343
4344 // Try to cast to a solid node
4345 SolidNode* const local_solid_node_pt =
4346 dynamic_cast<SolidNode*>(local_node_pt);
4347 // If we are successful
4348 if (local_solid_node_pt)
4349 {
4350 // Find the dimension of the node
4351 const unsigned n_dim = local_solid_node_pt->ndim();
4352 // Find number of positions
4353 const unsigned n_position_type =
4354 local_solid_node_pt->nposition_type();
4355
4356 for (unsigned k = 0; k < n_position_type; k++)
4357 {
4358 for (unsigned i = 0; i < n_dim; i++)
4359 {
4360 // If the generalised position is pinned,
4361 // set the value to zero
4362 if (local_solid_node_pt->position_is_pinned(k, i))
4363 {
4364 local_solid_node_pt->x_gen(k, i) = 0.0;
4365 }
4366 }
4367 }
4368 }
4369 }
4370
4371 // Now loop over the element's and zero the internal data
4372 const unsigned n_element = Mesh_pt->nelement();
4373 for (unsigned e = 0; e < n_element; e++)
4374 {
4375 GeneralisedElement* const local_element_pt = Mesh_pt->element_pt(e);
4376 const unsigned n_internal = local_element_pt->ninternal_data();
4377 for (unsigned i = 0; i < n_internal; i++)
4378 {
4379 Data* const local_data_pt = local_element_pt->internal_data_pt(i);
4380 const unsigned n_value = local_data_pt->nvalue();
4381 for (unsigned j = 0; j < n_value; j++)
4382 {
4383 // If the data value is pinned set the value to zero
4384 if (local_data_pt->is_pinned(j))
4385 {
4386 local_data_pt->set_value(j, 0.0);
4387 }
4388 }
4389 }
4390 } // End of loop over elements
4391 }
4392 else
4393 {
4394 // Alternatively loop over all sub meshes
4395 for (unsigned m = 0; m < n_sub_mesh; m++)
4396 {
4397 // Loop over the nodes in the element
4398 const unsigned n_node = Sub_mesh_pt[m]->nnode();
4399 for (unsigned n = 0; n < n_node; n++)
4400 {
4401 Node* const local_node_pt = Sub_mesh_pt[m]->node_pt(n);
4402 const unsigned n_value = local_node_pt->nvalue();
4403 for (unsigned j = 0; j < n_value; j++)
4404 {
4405 // If the data value is pinned set the value to zero
4406 if (local_node_pt->is_pinned(j))
4407 {
4408 local_node_pt->set_value(j, 0.0);
4409 }
4410 }
4411
4412 // Try to cast to a solid node
4413 SolidNode* const local_solid_node_pt =
4414 dynamic_cast<SolidNode*>(local_node_pt);
4415 // If we are successful
4416 if (local_solid_node_pt)
4417 {
4418 // Find the dimension of the node
4419 const unsigned n_dim = local_solid_node_pt->ndim();
4420 // Find number of positions
4421 const unsigned n_position_type =
4422 local_solid_node_pt->nposition_type();
4423
4424 for (unsigned k = 0; k < n_position_type; k++)
4425 {
4426 for (unsigned i = 0; i < n_dim; i++)
4427 {
4428 // If the generalised position is pinned,
4429 // set the value to zero
4430 if (local_solid_node_pt->position_is_pinned(k, i))
4431 {
4432 local_solid_node_pt->x_gen(k, i) = 0.0;
4433 }
4434 }
4435 }
4436 }
4437 }
4438
4439 // Now loop over the element's and zero the internal data
4440 const unsigned n_element = Sub_mesh_pt[m]->nelement();
4441 for (unsigned e = 0; e < n_element; e++)
4442 {
4443 GeneralisedElement* const local_element_pt =
4444 Sub_mesh_pt[m]->element_pt(e);
4445 const unsigned n_internal = local_element_pt->ninternal_data();
4446 for (unsigned i = 0; i < n_internal; i++)
4447 {
4448 Data* const local_data_pt = local_element_pt->internal_data_pt(i);
4449 const unsigned n_value = local_data_pt->nvalue();
4450 for (unsigned j = 0; j < n_value; j++)
4451 {
4452 // If the data value is pinned set the value to zero
4453 if (local_data_pt->is_pinned(j))
4454 {
4455 local_data_pt->set_value(j, 0.0);
4456 }
4457 }
4458 }
4459 } // End of loop over elements
4460 }
4461 }
4462 }
4463
4464
4465 //=====================================================================
4466 /// This is a (private) helper function that is used to assemble system
4467 /// matrices in compressed row or column format
4468 /// and compute residual vectors.
4469 /// The default action is to assemble the jacobian matrix and
4470 /// residuals for the Newton method. The action can be
4471 /// overloaded at an elemental level by changing the default
4472 /// behaviour of the function Element::get_all_vectors_and_matrices().
4473 /// column_or_row_index: Column [or row] index of given entry
4474 /// row_or_column_start: Index of first entry for given row [or column]
4475 /// value : Vector of nonzero entries
4476 /// residuals : Residual vector
4477 /// compressed_row_flag: Bool flag to indicate if storage format is
4478 /// compressed row [if false interpretation of
4479 /// arguments is as stated in square brackets].
4480 /// We provide four different assembly methods, each with different
4481 /// memory requirements/execution speeds. The method is set by
4482 /// the public flag Problem::Sparse_assembly_method.
4483 //=====================================================================
4485 Vector<int*>& column_or_row_index,
4486 Vector<int*>& row_or_column_start,
4487 Vector<double*>& value,
4488 Vector<unsigned>& nnz,
4489 Vector<double*>& residuals,
4490 bool compressed_row_flag)
4491 {
4492 // Choose the actual method
4493 switch (Sparse_assembly_method)
4494 {
4496
4498 column_or_row_index,
4499 row_or_column_start,
4500 value,
4501 nnz,
4502 residuals,
4503 compressed_row_flag);
4504
4505 break;
4506
4508
4510 column_or_row_index,
4511 row_or_column_start,
4512 value,
4513 nnz,
4514 residuals,
4515 compressed_row_flag);
4516
4517 break;
4518
4520
4522 row_or_column_start,
4523 value,
4524 nnz,
4525 residuals,
4526 compressed_row_flag);
4527
4528 break;
4529
4531
4533 column_or_row_index,
4534 row_or_column_start,
4535 value,
4536 nnz,
4537 residuals,
4538 compressed_row_flag);
4539
4540 break;
4541
4543
4545 column_or_row_index,
4546 row_or_column_start,
4547 value,
4548 nnz,
4549 residuals,
4550 compressed_row_flag);
4551
4552 break;
4553
4554 default:
4555
4556 std::ostringstream error_stream;
4557 error_stream
4558 << "Error: Incorrect value for Problem::Sparse_assembly_method"
4559 << Sparse_assembly_method << std::endl
4560 << "It should be one of the enumeration Problem::Assembly_method"
4561 << std::endl;
4562 throw OomphLibError(
4563 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4564 }
4565 }
4566
4567
4568 //=====================================================================
4569 /// This is a (private) helper function that is used to assemble system
4570 /// matrices in compressed row or column format
4571 /// and compute residual vectors, using maps
4572 /// The default action is to assemble the jacobian matrix and
4573 /// residuals for the Newton method. The action can be
4574 /// overloaded at an elemental level by chaging the default
4575 /// behaviour of the function Element::get_all_vectors_and_matrices().
4576 /// column_or_row_index: Column [or row] index of given entry
4577 /// row_or_column_start: Index of first entry for given row [or column]
4578 /// value : Vector of nonzero entries
4579 /// residuals : Residual vector
4580 /// compressed_row_flag: Bool flag to indicate if storage format is
4581 /// compressed row [if false interpretation of
4582 /// arguments is as stated in square brackets].
4583 //=====================================================================
4585 Vector<int*>& column_or_row_index,
4586 Vector<int*>& row_or_column_start,
4587 Vector<double*>& value,
4588 Vector<unsigned>& nnz,
4589 Vector<double*>& residuals,
4590 bool compressed_row_flag)
4591 {
4592 // Total number of elements
4593 const unsigned long n_elements = mesh_pt()->nelement();
4594
4595 // Default range of elements for distributed problems
4596 unsigned long el_lo = 0;
4597 unsigned long el_hi = n_elements - 1;
4598
4599#ifdef OOMPH_HAS_MPI
4600 // Otherwise just loop over a fraction of the elements
4601 // (This will either have been initialised in
4602 // Problem::set_default_first_and_last_element_for_assembly() or
4603 // will have been re-assigned during a previous assembly loop
4604 // Note that following the re-assignment only the entries
4605 // for the current processor are relevant.
4607 {
4608 el_lo = First_el_for_assembly[Communicator_pt->my_rank()];
4609 el_hi = Last_el_plus_one_for_assembly[Communicator_pt->my_rank()] - 1;
4610 }
4611#endif
4612
4613 // number of dofs
4614 unsigned ndof = this->ndof();
4615
4616 // Find the number of vectors to be assembled
4617 const unsigned n_vector = residuals.size();
4618
4619 // Find the number of matrices to be assembled
4620 const unsigned n_matrix = column_or_row_index.size();
4621
4622 // Locally cache pointer to assembly handler
4624
4625#ifdef OOMPH_HAS_MPI
4626 bool doing_residuals = false;
4627 if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt) != 0)
4628 {
4629 doing_residuals = true;
4630 }
4631#endif
4632
4633// Error check dimensions
4634#ifdef PARANOID
4635 if (row_or_column_start.size() != n_matrix)
4636 {
4637 std::ostringstream error_stream;
4638 error_stream << "Error: " << std::endl
4639 << "row_or_column_start.size() "
4640 << row_or_column_start.size() << " does not equal "
4641 << "column_or_row_index.size() "
4642 << column_or_row_index.size() << std::endl;
4643 throw OomphLibError(
4644 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4645 }
4646
4647 if (value.size() != n_matrix)
4648 {
4649 std::ostringstream error_stream;
4650 error_stream
4651 << "Error in Problem::sparse_assemble_row_or_column_compressed "
4652 << std::endl
4653 << "value.size() " << value.size() << " does not equal "
4654 << "column_or_row_index.size() " << column_or_row_index.size()
4655 << std::endl
4656 << std::endl
4657 << std::endl;
4658 throw OomphLibError(
4659 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4660 }
4661#endif
4662
4663
4664 // The idea behind this sparse assembly routine is to use a vector of
4665 // maps for the entries in each row or column of the complete matrix.
4666 // The key for each map is the global row or column number and
4667 // the default comparison operator for integers means that each map
4668 // is ordered by the global row or column number. Thus, we need not
4669 // sort the maps, that happens at each insertion of a new entry. The
4670 // price we pay is that for large maps, inseration is not a
4671 // cheap operation. Hash maps can be used to increase the speed, but then
4672 // the ordering is lost and we would have to sort anyway. The solution if
4673 // speed is required is to use lists, see below.
4674
4675
4676 // Set up a vector of vectors of maps of entries of each matrix,
4677 // indexed by either the column or row. The entries of the vector for
4678 // each matrix correspond to all the rows or columns of that matrix.
4679 // The use of the map storage
4680 // scheme, with its implicit ordering on the first index, gives
4681 // a sparse ordered list of the entries in the given row or column.
4682 Vector<Vector<std::map<unsigned, double>>> matrix_data_map(n_matrix);
4683 // Loop over the number of matrices being assembled and resize
4684 // each vector of maps to the number of rows or columns of the matrix
4685 for (unsigned m = 0; m < n_matrix; m++)
4686 {
4687 matrix_data_map[m].resize(ndof);
4688 }
4689
4690 // Resize the residuals vectors
4691 for (unsigned v = 0; v < n_vector; v++)
4692 {
4693 residuals[v] = new double[ndof];
4694 for (unsigned i = 0; i < ndof; i++)
4695 {
4696 residuals[v][i] = 0;
4697 }
4698 }
4699
4700
4701#ifdef OOMPH_HAS_MPI
4702
4703
4704 // Storage for assembly time for elements
4705 double t_assemble_start = 0.0;
4706
4707 // Storage for assembly times
4708 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
4709 {
4710 Elemental_assembly_time.resize(n_elements);
4711 }
4712
4713#endif
4714
4715 //----------------Assemble and populate the maps-------------------------
4716 {
4717 // Allocate local storage for the element's contribution to the
4718 // residuals vectors and system matrices of the size of the maximum
4719 // number of dofs in any element.
4720 // This means that the storage is only allocated (and deleted) once
4721 Vector<Vector<double>> el_residuals(n_vector);
4722 Vector<DenseMatrix<double>> el_jacobian(n_matrix);
4723
4724 // Loop over the elements for this processor
4725 for (unsigned long e = el_lo; e <= el_hi; e++)
4726 {
4727#ifdef OOMPH_HAS_MPI
4728 // Time it?
4729 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
4730 {
4731 t_assemble_start = TimingHelpers::timer();
4732 }
4733#endif
4734
4735 // Get the pointer to the element
4736 GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
4737
4738#ifdef OOMPH_HAS_MPI
4739 // Ignore halo elements
4740 if (!elem_pt->is_halo())
4741 {
4742#endif
4743
4744 // Find number of degrees of freedom in the element
4745 const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
4746
4747 // Resize the storage for elemental jacobian and residuals
4748 for (unsigned v = 0; v < n_vector; v++)
4749 {
4750 el_residuals[v].resize(nvar);
4751 }
4752 for (unsigned m = 0; m < n_matrix; m++)
4753 {
4754 el_jacobian[m].resize(nvar);
4755 }
4756
4757 // Now get the residuals and jacobian for the element
4759 elem_pt, el_residuals, el_jacobian);
4760
4761 //---------------Insert the values into the maps--------------
4762
4763 // Loop over the first index of local variables
4764 for (unsigned i = 0; i < nvar; i++)
4765 {
4766 // Get the local equation number
4767 unsigned eqn_number = assembly_handler_pt->eqn_number(elem_pt, i);
4768
4769 // Add the contribution to the residuals
4770 for (unsigned v = 0; v < n_vector; v++)
4771 {
4772 // Fill in each residuals vector
4773 residuals[v][eqn_number] += el_residuals[v][i];
4774 }
4775
4776 // Now loop over the other index
4777 for (unsigned j = 0; j < nvar; j++)
4778 {
4779 // Get the number of the unknown
4780 unsigned unknown = assembly_handler_pt->eqn_number(elem_pt, j);
4781
4782 // Loop over the matrices
4783 for (unsigned m = 0; m < n_matrix; m++)
4784 {
4785 // Get the value of the matrix at this point
4786 double value = el_jacobian[m](i, j);
4787 // Only bother to add to the map if it's non-zero
4788 if (std::fabs(value) > Numerical_zero_for_sparse_assembly)
4789 {
4790 // If it's compressed row storage, then our vector of maps
4791 // is indexed by row (equation number)
4792 if (compressed_row_flag)
4793 {
4794 // Add the data into the map using the unknown as the map
4795 // key
4796 matrix_data_map[m][eqn_number][unknown] += value;
4797 }
4798 // Otherwise it's compressed column storage and our vector is
4799 // indexed by column (the unknown)
4800 else
4801 {
4802 // Add the data into the map using the eqn_numbe as the map
4803 // key
4804 matrix_data_map[m][unknown][eqn_number] += value;
4805 }
4806 }
4807 } // End of loop over matrices
4808 }
4809 }
4810
4811#ifdef OOMPH_HAS_MPI
4812 } // endif halo element
4813#endif
4814
4815
4816#ifdef OOMPH_HAS_MPI
4817 // Time it?
4818 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
4819 {
4821 TimingHelpers::timer() - t_assemble_start;
4822 }
4823#endif
4824
4825 } // End of loop over the elements
4826
4827 } // End of map assembly
4828
4829
4830#ifdef OOMPH_HAS_MPI
4831
4832 // Postprocess timing information and re-allocate distribution of
4833 // elements during subsequent assemblies.
4834 if ((!doing_residuals) && (!Problem_has_been_distributed) &&
4836 {
4838 }
4839
4840 // We have determined load balancing for current setup.
4841 // This can remain the same until assign_eqn_numbers() is called
4842 // again -- the flag is re-set to true there.
4843 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
4844 {
4846 }
4847
4848#endif
4849
4850
4851 //-----------Finally we need to convert the beautiful map storage scheme
4852 //------------------------to the containers required by SuperLU
4853
4854 // Loop over the number of matrices
4855 for (unsigned m = 0; m < n_matrix; m++)
4856 {
4857 // Set the number of rows or columns
4858 row_or_column_start[m] = new int[ndof + 1];
4859 // Counter for the total number of entries in the storage scheme
4860 unsigned long entry_count = 0;
4861 row_or_column_start[m][0] = entry_count;
4862
4863 // first we compute the number of non-zeros
4864 nnz[m] = 0;
4865 for (unsigned long i_global = 0; i_global < ndof; i_global++)
4866 {
4867 nnz[m] += matrix_data_map[m][i_global].size();
4868 }
4869
4870 // and then resize the storage
4871 column_or_row_index[m] = new int[nnz[m]];
4872 value[m] = new double[nnz[m]];
4873
4874 // Now we merely loop over the number of rows or columns
4875 for (unsigned long i_global = 0; i_global < ndof; i_global++)
4876 {
4877 // Start index for the present row
4878 row_or_column_start[m][i_global] = entry_count;
4879 // If there are no entries in the map then skip the rest of the loop
4880 if (matrix_data_map[m][i_global].empty())
4881 {
4882 continue;
4883 }
4884
4885 // Loop over all the entries in the map corresponding to the given
4886 // row or column. It will be ordered
4887
4888 for (std::map<unsigned, double>::iterator it =
4889 matrix_data_map[m][i_global].begin();
4890 it != matrix_data_map[m][i_global].end();
4891 ++it)
4892 {
4893 // The first value is the column or row index
4894 column_or_row_index[m][entry_count] = it->first;
4895 // The second value is the actual data value
4896 value[m][entry_count] = it->second;
4897 // Increase the value of the counter
4898 entry_count++;
4899 }
4900 }
4901
4902 // Final entry in the row/column start vector
4903 row_or_column_start[m][ndof] = entry_count;
4904 } // End of the loop over the matrices
4905
4907 {
4908 oomph_info << "Pausing at end of sparse assembly." << std::endl;
4909 pause("Check memory usage now.");
4910 }
4911 }
4912
4913
4914 //=====================================================================
4915 /// This is a (private) helper function that is used to assemble system
4916 /// matrices in compressed row or column format
4917 /// and compute residual vectors using lists
4918 /// The default action is to assemble the jacobian matrix and
4919 /// residuals for the Newton method. The action can be
4920 /// overloaded at an elemental level by chaging the default
4921 /// behaviour of the function Element::get_all_vectors_and_matrices().
4922 /// column_or_row_index: Column [or row] index of given entry
4923 /// row_or_column_start: Index of first entry for given row [or column]
4924 /// value : Vector of nonzero entries
4925 /// residuals : Residual vector
4926 /// compressed_row_flag: Bool flag to indicate if storage format is
4927 /// compressed row [if false interpretation of
4928 /// arguments is as stated in square brackets].
4929 //=====================================================================
4931 Vector<int*>& column_or_row_index,
4932 Vector<int*>& row_or_column_start,
4933 Vector<double*>& value,
4934 Vector<unsigned>& nnz,
4935 Vector<double*>& residuals,
4936 bool compressed_row_flag)
4937 {
4938 // Total number of elements
4939 const unsigned long n_elements = mesh_pt()->nelement();
4940
4941 // Default range of elements for distributed problems
4942 unsigned long el_lo = 0;
4943 unsigned long el_hi = n_elements - 1;
4944
4945#ifdef OOMPH_HAS_MPI
4946 // Otherwise just loop over a fraction of the elements
4947 // (This will either have been initialised in
4948 // Problem::set_default_first_and_last_element_for_assembly() or
4949 // will have been re-assigned during a previous assembly loop
4950 // Note that following the re-assignment only the entries
4951 // for the current processor are relevant.
4953 {
4954 el_lo = First_el_for_assembly[Communicator_pt->my_rank()];
4955 el_hi = Last_el_plus_one_for_assembly[Communicator_pt->my_rank()] - 1;
4956 }
4957#endif
4958
4959 // number of dofs
4960 unsigned ndof = this->ndof();
4961
4962 // Find the number of vectors to be assembled
4963 const unsigned n_vector = residuals.size();
4964
4965 // Find the number of matrices to be assembled
4966 const unsigned n_matrix = column_or_row_index.size();
4967
4968 // Locally cache pointer to assembly handler
4970
4971#ifdef OOMPH_HAS_MPI
4972 bool doing_residuals = false;
4973 if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt) != 0)
4974 {
4975 doing_residuals = true;
4976 }
4977#endif
4978
4979// Error check dimensions
4980#ifdef PARANOID
4981 if (row_or_column_start.size() != n_matrix)
4982 {
4983 std::ostringstream error_stream;
4984 error_stream << "Error: " << std::endl
4985 << "row_or_column_start.size() "
4986 << row_or_column_start.size() << " does not equal "
4987 << "column_or_row_index.size() "
4988 << column_or_row_index.size() << std::endl;
4989 throw OomphLibError(
4990 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4991 }
4992
4993 if (value.size() != n_matrix)
4994 {
4995 std::ostringstream error_stream;
4996 error_stream
4997 << "Error in Problem::sparse_assemble_row_or_column_compressed "
4998 << std::endl
4999 << "value.size() " << value.size() << " does not equal "
5000 << "column_or_row_index.size() " << column_or_row_index.size()
5001 << std::endl
5002 << std::endl
5003 << std::endl;
5004 throw OomphLibError(
5005 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5006 }
5007#endif
5008
5009 // The idea behind this sparse assembly routine is to use a vector of
5010 // lists for the entries in each row or column of the complete matrix.
5011 // The lists contain pairs of entries (global row/column number, value).
5012 // All non-zero contributions from each element are added to the lists.
5013 // We then sort each list by global row/column number and then combine
5014 // the entries corresponding to each row/column before adding to the
5015 // vectors column_or_row_index and value.
5016
5017 // Note the trade off for "fast assembly" is that we will require
5018 // more memory during the assembly phase. Then again, if we can
5019 // only just assemble the sparse matrix, we're in real trouble.
5020
5021 // Set up a vector of lists of paired entries of
5022 //(row/column index, jacobian matrix entry).
5023 // The entries of the vector correspond to all the rows or columns.
5024 // The use of the list storage scheme, should give fast insertion
5025 // and fast sorts later.
5027 n_matrix);
5028 // Loop over the number of matrices and resize
5029 for (unsigned m = 0; m < n_matrix; m++)
5030 {
5031 matrix_data_list[m].resize(ndof);
5032 }
5033
5034 // Resize the residuals vectors
5035 for (unsigned v = 0; v < n_vector; v++)
5036 {
5037 residuals[v] = new double[ndof];
5038 for (unsigned i = 0; i < ndof; i++)
5039 {
5040 residuals[v][i] = 0;
5041 }
5042 }
5043
5044#ifdef OOMPH_HAS_MPI
5045
5046
5047 // Storage for assembly time for elements
5048 double t_assemble_start = 0.0;
5049
5050 // Storage for assembly times
5051 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5052 {
5053 Elemental_assembly_time.resize(n_elements);
5054 }
5055
5056#endif
5057
5058 //------------Assemble and populate the lists-----------------------
5059 {
5060 // Allocate local storage for the element's contribution to the
5061 // residuals vectors and system matrices of the size of the maximum
5062 // number of dofs in any element.
5063 // This means that the stored is only allocated (and deleted) once
5064 Vector<Vector<double>> el_residuals(n_vector);
5065 Vector<DenseMatrix<double>> el_jacobian(n_matrix);
5066
5067
5068 // Pointer to a single list to be used during the assembly
5069 std::list<std::pair<unsigned, double>>* list_pt;
5070
5071 // Loop over the all elements
5072 for (unsigned long e = el_lo; e <= el_hi; e++)
5073 {
5074#ifdef OOMPH_HAS_MPI
5075 // Time it?
5076 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5077 {
5078 t_assemble_start = TimingHelpers::timer();
5079 }
5080#endif
5081
5082 // Get the pointer to the element
5083 GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
5084
5085#ifdef OOMPH_HAS_MPI
5086 // Ignore halo elements
5087 if (!elem_pt->is_halo())
5088 {
5089#endif
5090
5091 // Find number of degrees of freedom in the element
5092 const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
5093
5094 // Resize the storage for the elemental jacobian and residuals
5095 for (unsigned v = 0; v < n_vector; v++)
5096 {
5097 el_residuals[v].resize(nvar);
5098 }
5099 for (unsigned m = 0; m < n_matrix; m++)
5100 {
5101 el_jacobian[m].resize(nvar);
5102 }
5103
5104 // Now get the residuals and jacobian for the element
5106 elem_pt, el_residuals, el_jacobian);
5107
5108 //---------------- Insert the values into the lists -----------
5109
5110 // Loop over the first index of local variables
5111 for (unsigned i = 0; i < nvar; i++)
5112 {
5113 // Get the local equation number
5114 unsigned eqn_number = assembly_handler_pt->eqn_number(elem_pt, i);
5115
5116 // Add the contribution to the residuals
5117 for (unsigned v = 0; v < n_vector; v++)
5118 {
5119 // Fill in the residuals vector
5120 residuals[v][eqn_number] += el_residuals[v][i];
5121 }
5122
5123 // Now loop over the other index
5124 for (unsigned j = 0; j < nvar; j++)
5125 {
5126 // Get the number of the unknown
5127 unsigned unknown = assembly_handler_pt->eqn_number(elem_pt, j);
5128
5129 // Loop over the matrices
5130 for (unsigned m = 0; m < n_matrix; m++)
5131 {
5132 // Get the value of the matrix at this point
5133 double value = el_jacobian[m](i, j);
5134 // Only add to theif it's non-zero
5135 if (std::fabs(value) > Numerical_zero_for_sparse_assembly)
5136 {
5137 // If it's compressed row storage, then our vector is indexed
5138 // by row (the equation number)
5139 if (compressed_row_flag)
5140 {
5141 // Find the list that corresponds to the desired row
5142 list_pt = &matrix_data_list[m][eqn_number];
5143 // Insert the data into the list, the first entry
5144 // in the pair is the unknown (column index),
5145 // the second is the value itself.
5146 list_pt->insert(list_pt->end(),
5147 std::make_pair(unknown, value));
5148 }
5149 // Otherwise it's compressed column storage, and our
5150 // vector is indexed by column (the unknown)
5151 else
5152 {
5153 // Find the list that corresponds to the desired column
5154 list_pt = &matrix_data_list[m][unknown];
5155 // Insert the data into the list, the first entry
5156 // in the pair is the equation number (row index),
5157 // the second is the value itself.
5158 list_pt->insert(list_pt->end(),
5159 std::make_pair(eqn_number, value));
5160 }
5161 }
5162 }
5163 }
5164 }
5165
5166#ifdef OOMPH_HAS_MPI
5167 } // endif halo element
5168#endif
5169
5170
5171#ifdef OOMPH_HAS_MPI
5172 // Time it?
5173 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5174 {
5176 TimingHelpers::timer() - t_assemble_start;
5177 }
5178#endif
5179
5180 } // End of loop over the elements
5181
5182 } // list_pt goes out of scope
5183
5184
5185#ifdef OOMPH_HAS_MPI
5186
5187 // Postprocess timing information and re-allocate distribution of
5188 // elements during subsequent assemblies.
5189 if ((!doing_residuals) && (!Problem_has_been_distributed) &&
5191 {
5193 }
5194
5195 // We have determined load balancing for current setup.
5196 // This can remain the same until assign_eqn_numbers() is called
5197 // again -- the flag is re-set to true there.
5198 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5199 {
5201 }
5202
5203#endif
5204
5205
5206 //----Finally we need to convert the beautiful list storage scheme---
5207 //----------to the containers required by SuperLU--------------------
5208
5209 // Loop over the number of matrices
5210 for (unsigned m = 0; m < n_matrix; m++)
5211 {
5212 // Set the number of rows or columns
5213 row_or_column_start[m] = new int[ndof + 1];
5214 // Counter for the total number of entries in the storage scheme
5215 unsigned long entry_count = 0;
5216 // The first entry is 0
5217 row_or_column_start[m][0] = entry_count;
5218
5219 // first we compute the number of non-zeros
5220 nnz[m] = 0;
5221 for (unsigned long i_global = 0; i_global < ndof; i_global++)
5222 {
5223 nnz[m] += matrix_data_list[m][i_global].size();
5224 }
5225
5226 // and then resize the storage
5227 column_or_row_index[m] = new int[nnz[m]];
5228 value[m] = new double[nnz[m]];
5229
5230 // Now we merely loop over the number of rows or columns
5231 for (unsigned long i_global = 0; i_global < ndof; i_global++)
5232 {
5233 // Start index for the present row is the number of entries so far
5234 row_or_column_start[m][i_global] = entry_count;
5235 // If there are no entries in the list then skip the loop
5236 if (matrix_data_list[m][i_global].empty())
5237 {
5238 continue;
5239 }
5240
5241 // Sort the list corresponding to this row or column by the
5242 // column or row index (first entry in the pair).
5243 // This might be inefficient, but we only have to do the sort ONCE
5244 // for each list. This is faster than using a map storage scheme, where
5245 // we are sorting for every insertion (although the map structure
5246 // is cleaner and more memory efficient)
5247 matrix_data_list[m][i_global].sort();
5248
5249 // Set up an iterator for start of the list
5250 std::list<std::pair<unsigned, double>>::iterator it =
5251 matrix_data_list[m][i_global].begin();
5252
5253 // Get the first row or column index in the list...
5254 unsigned current_index = it->first;
5255 //...and the corresponding value
5256 double current_value = it->second;
5257
5258 // Loop over all the entries in the sorted list
5259 // Increase the iterator so that we start at the second entry
5260 for (++it; it != matrix_data_list[m][i_global].end(); ++it)
5261 {
5262 // If the index has not changed, then we must add the contribution
5263 // of the present entry to the value.
5264 // Additionally check that the entry is non-zero
5265 if ((it->first == current_index) &&
5266 (std::fabs(it->second) > Numerical_zero_for_sparse_assembly))
5267 {
5268 current_value += it->second;
5269 }
5270 // Otherwise, we have added all the contributions to the index
5271 // to current_value, so add it to the SuperLU data structure
5272 else
5273 {
5274 // Add the row or column index to the vector
5275 column_or_row_index[m][entry_count] = current_index;
5276 // Add the actual value to the vector
5277 value[m][entry_count] = current_value;
5278 // Increase the counter for the number of entries in each vector
5279 entry_count++;
5280
5281 // Set the index and value to be those of the current entry in the
5282 // list
5283 current_index = it->first;
5284 current_value = it->second;
5285 }
5286 } // End of loop over all list entries for this global row or column
5287
5288 // There are TWO special cases to consider.
5289 // If there is only one equation number in the list, then it
5290 // will NOT have been added. We test this case by comparing the
5291 // number of entries with those stored in row_or_column_start[i_global]
5292 // Otherwise
5293 // If the final entry in the list has the same index as the penultimate
5294 // entry, then it will NOT have been added to the SuperLU storage scheme
5295 // Check this by comparing the current_index with the final index
5296 // stored in the SuperLU scheme. If they are not the same, then
5297 // add the current_index and value.
5298
5299 // If single equation number in list
5300 if ((static_cast<int>(entry_count) == row_or_column_start[m][i_global])
5301 // If we have a single equation number, this will not be evaluated.
5302 // If we don't then we do the test to check that the final
5303 // entry is added
5304 || (static_cast<int>(current_index) !=
5305 column_or_row_index[m][entry_count - 1]))
5306 {
5307 // Add the row or column index to the vector
5308 column_or_row_index[m][entry_count] = current_index;
5309 // Add the actual value to the vector
5310 value[m][entry_count] = current_value;
5311 // Increase the counter for the number of entries in each vector
5312 entry_count++;
5313 }
5314
5315 } // End of loop over the rows or columns of the entire matrix
5316
5317 // Final entry in the row/column start vector
5318 row_or_column_start[m][ndof] = entry_count;
5319 } // End of loop over matrices
5320
5322 {
5323 oomph_info << "Pausing at end of sparse assembly." << std::endl;
5324 pause("Check memory usage now.");
5325 }
5326 }
5327
5328
5329 //=====================================================================
5330 /// This is a (private) helper function that is used to assemble system
5331 /// matrices in compressed row or column format
5332 /// and compute residual vectors using vectors of pairs
5333 /// The default action is to assemble the jacobian matrix and
5334 /// residuals for the Newton method. The action can be
5335 /// overloaded at an elemental level by chaging the default
5336 /// behaviour of the function Element::get_all_vectors_and_matrices().
5337 /// column_or_row_index: Column [or row] index of given entry
5338 /// row_or_column_start: Index of first entry for given row [or column]
5339 /// value : Vector of nonzero entries
5340 /// residuals : Residual vector
5341 /// compressed_row_flag: Bool flag to indicate if storage format is
5342 /// compressed row [if false interpretation of
5343 /// arguments is as stated in square brackets].
5344 //=====================================================================
5346 Vector<int*>& column_or_row_index,
5347 Vector<int*>& row_or_column_start,
5348 Vector<double*>& value,
5349 Vector<unsigned>& nnz,
5350 Vector<double*>& residuals,
5351 bool compressed_row_flag)
5352 {
5353 // Total number of elements
5354 const unsigned long n_elements = mesh_pt()->nelement();
5355
5356 // Default range of elements for distributed problems
5357 unsigned long el_lo = 0;
5358 unsigned long el_hi = n_elements - 1;
5359
5360#ifdef OOMPH_HAS_MPI
5361 // Otherwise just loop over a fraction of the elements
5362 // (This will either have been initialised in
5363 // Problem::set_default_first_and_last_element_for_assembly() or
5364 // will have been re-assigned during a previous assembly loop
5365 // Note that following the re-assignment only the entries
5366 // for the current processor are relevant.
5368 {
5369 el_lo = First_el_for_assembly[Communicator_pt->my_rank()];
5370 el_hi = Last_el_plus_one_for_assembly[Communicator_pt->my_rank()] - 1;
5371 }
5372#endif
5373
5374 // number of local eqns
5375 unsigned ndof = this->ndof();
5376
5377 // Find the number of vectors to be assembled
5378 const unsigned n_vector = residuals.size();
5379
5380 // Find the number of matrices to be assembled
5381 const unsigned n_matrix = column_or_row_index.size();
5382
5383 // Locally cache pointer to assembly handler
5385
5386#ifdef OOMPH_HAS_MPI
5387 bool doing_residuals = false;
5388 if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt) != 0)
5389 {
5390 doing_residuals = true;
5391 }
5392#endif
5393
5394// Error check dimensions
5395#ifdef PARANOID
5396 if (row_or_column_start.size() != n_matrix)
5397 {
5398 std::ostringstream error_stream;
5399 error_stream << "Error: " << std::endl
5400 << "row_or_column_start.size() "
5401 << row_or_column_start.size() << " does not equal "
5402 << "column_or_row_index.size() "
5403 << column_or_row_index.size() << std::endl;
5404 throw OomphLibError(
5405 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5406 }
5407
5408 if (value.size() != n_matrix)
5409 {
5410 std::ostringstream error_stream;
5411 error_stream << "Error: " << std::endl
5412 << "value.size() " << value.size() << " does not equal "
5413 << "column_or_row_index.size() "
5414 << column_or_row_index.size() << std::endl
5415 << std::endl
5416 << std::endl;
5417 throw OomphLibError(
5418 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5419 }
5420#endif
5421
5422
5423 // The idea behind this sparse assembly routine is to use a Vector of
5424 // Vectors of pairs for each complete matrix.
5425 // Each inner Vector stores pairs and holds the row (or column) index
5426 // and the value of the matrix entry.
5427
5428 // Set up Vector of Vectors to store the entries of each matrix,
5429 // indexed by either the column or row.
5430 Vector<Vector<Vector<std::pair<unsigned, double>>>> matrix_data(n_matrix);
5431
5432 // Loop over the number of matrices being assembled and resize
5433 // each Vector of Vectors to the number of rows or columns of the matrix
5434 for (unsigned m = 0; m < n_matrix; m++)
5435 {
5436 matrix_data[m].resize(ndof);
5437 }
5438
5439 // Resize the residuals vectors
5440 for (unsigned v = 0; v < n_vector; v++)
5441 {
5442 residuals[v] = new double[ndof];
5443 for (unsigned i = 0; i < ndof; i++)
5444 {
5445 residuals[v][i] = 0;
5446 }
5447 }
5448
5449#ifdef OOMPH_HAS_MPI
5450
5451 // Storage for assembly time for elements
5452 double t_assemble_start = 0.0;
5453
5454 // Storage for assembly times
5455 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5456 {
5457 Elemental_assembly_time.resize(n_elements);
5458 }
5459
5460#endif
5461
5462 //----------------Assemble and populate the vector storage scheme--------
5463 {
5464 // Allocate local storage for the element's contribution to the
5465 // residuals vectors and system matrices of the size of the maximum
5466 // number of dofs in any element
5467 // This means that the storage is only allocated (and deleted) once
5468 Vector<Vector<double>> el_residuals(n_vector);
5469 Vector<DenseMatrix<double>> el_jacobian(n_matrix);
5470
5471 // Loop over the elements
5472 for (unsigned long e = el_lo; e <= el_hi; e++)
5473 {
5474#ifdef OOMPH_HAS_MPI
5475 // Time it?
5476 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5477 {
5478 t_assemble_start = TimingHelpers::timer();
5479 }
5480#endif
5481
5482 // Get the pointer to the element
5483 GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
5484
5485#ifdef OOMPH_HAS_MPI
5486 // Ignore halo elements
5487 if (!elem_pt->is_halo())
5488 {
5489#endif
5490
5491 // Find number of degrees of freedom in the element
5492 const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
5493
5494 // Resize the storage for elemental jacobian and residuals
5495 for (unsigned v = 0; v < n_vector; v++)
5496 {
5497 el_residuals[v].resize(nvar);
5498 }
5499 for (unsigned m = 0; m < n_matrix; m++)
5500 {
5501 el_jacobian[m].resize(nvar);
5502 }
5503
5504 // Now get the residuals and jacobian for the element
5506 elem_pt, el_residuals, el_jacobian);
5507
5508 //---------------Insert the values into the vectors--------------
5509
5510 // Loop over the first index of local variables
5511 for (unsigned i = 0; i < nvar; i++)
5512 {
5513 // Get the local equation number
5514 unsigned eqn_number = assembly_handler_pt->eqn_number(elem_pt, i);
5515
5516 // Add the contribution to the residuals
5517 for (unsigned v = 0; v < n_vector; v++)
5518 {
5519 // Fill in each residuals vector
5520 residuals[v][eqn_number] += el_residuals[v][i];
5521 }
5522
5523 // Now loop over the other index
5524 for (unsigned j = 0; j < nvar; j++)
5525 {
5526 // Get the number of the unknown
5527 unsigned unknown = assembly_handler_pt->eqn_number(elem_pt, j);
5528
5529 // Loop over the matrices
5530 // If it's compressed row storage, then our vector of maps
5531 // is indexed by row (equation number)
5532 for (unsigned m = 0; m < n_matrix; m++)
5533 {
5534 // Get the value of the matrix at this point
5535 double value = el_jacobian[m](i, j);
5536 // Only bother to add to the vector if it's non-zero
5537 if (std::fabs(value) > Numerical_zero_for_sparse_assembly)
5538 {
5539 // If it's compressed row storage, then our vector of maps
5540 // is indexed by row (equation number)
5541 if (compressed_row_flag)
5542 {
5543 // Find the correct position and add the data into the
5544 // vectors
5545 const unsigned size = matrix_data[m][eqn_number].size();
5546 for (unsigned k = 0; k <= size; k++)
5547 {
5548 if (k == size)
5549 {
5550 matrix_data[m][eqn_number].push_back(
5551 std::make_pair(unknown, value));
5552 break;
5553 }
5554 else if (matrix_data[m][eqn_number][k].first == unknown)
5555 {
5556 matrix_data[m][eqn_number][k].second += value;
5557 break;
5558 }
5559 }
5560 }
5561 // Otherwise it's compressed column storage and our vector is
5562 // indexed by column (the unknown)
5563 else
5564 {
5565 // Add the data into the vectors in the correct position
5566 const unsigned size = matrix_data[m][unknown].size();
5567 for (unsigned k = 0; k <= size; k++)
5568 {
5569 if (k == size)
5570 {
5571 matrix_data[m][unknown].push_back(
5572 std::make_pair(eqn_number, value));
5573 break;
5574 }
5575 else if (matrix_data[m][unknown][k].first == eqn_number)
5576 {
5577 matrix_data[m][unknown][k].second += value;
5578 break;
5579 }
5580 }
5581 }
5582 }
5583 } // End of loop over matrices
5584 }
5585 }
5586
5587#ifdef OOMPH_HAS_MPI
5588 } // endif halo element
5589#endif
5590
5591
5592#ifdef OOMPH_HAS_MPI
5593 // Time it?
5594 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5595 {
5597 TimingHelpers::timer() - t_assemble_start;
5598 }
5599#endif
5600
5601 } // End of loop over the elements
5602
5603
5604 } // End of vector assembly
5605
5606
5607#ifdef OOMPH_HAS_MPI
5608
5609 // Postprocess timing information and re-allocate distribution of
5610 // elements during subsequent assemblies.
5611 if ((!doing_residuals) && (!Problem_has_been_distributed) &&
5613 {
5615 }
5616
5617 // We have determined load balancing for current setup.
5618 // This can remain the same until assign_eqn_numbers() is called
5619 // again -- the flag is re-set to true there.
5620 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5621 {
5623 }
5624
5625#endif
5626
5627
5628 //-----------Finally we need to convert this vector storage scheme
5629 //------------------------to the containers required by SuperLU
5630
5631 // Loop over the number of matrices
5632 for (unsigned m = 0; m < n_matrix; m++)
5633 {
5634 // Set the number of rows or columns
5635 row_or_column_start[m] = new int[ndof + 1];
5636
5637 // fill row_or_column_start and find the number of entries
5638 row_or_column_start[m][0] = 0;
5639 for (unsigned long i = 0; i < ndof; i++)
5640 {
5641 row_or_column_start[m][i + 1] =
5642 row_or_column_start[m][i] + matrix_data[m][i].size();
5643 }
5644 const unsigned entries = row_or_column_start[m][ndof];
5645
5646 // resize vectors
5647 column_or_row_index[m] = new int[entries];
5648 value[m] = new double[entries];
5649 nnz[m] = entries;
5650
5651 // Now we merely loop over the number of rows or columns
5652 for (unsigned long i_global = 0; i_global < ndof; i_global++)
5653 {
5654 // If there are no entries in the vector then skip the rest of the loop
5655 if (matrix_data[m][i_global].empty())
5656 {
5657 continue;
5658 }
5659
5660 // Loop over all the entries in the vectors corresponding to the given
5661 // row or column. It will NOT be ordered
5662 unsigned p = 0;
5663 for (int j = row_or_column_start[m][i_global];
5664 j < row_or_column_start[m][i_global + 1];
5665 j++)
5666 {
5667 column_or_row_index[m][j] = matrix_data[m][i_global][p].first;
5668 value[m][j] = matrix_data[m][i_global][p].second;
5669 ++p;
5670 }
5671 }
5672 } // End of the loop over the matrices
5673
5675 {
5676 oomph_info << "Pausing at end of sparse assembly." << std::endl;
5677 pause("Check memory usage now.");
5678 }
5679 }
5680
5681
5682 //=====================================================================
5683 /// This is a (private) helper function that is used to assemble system
5684 /// matrices in compressed row or column format
5685 /// and compute residual vectors using two vectors.
5686 /// The default action is to assemble the jacobian matrix and
5687 /// residuals for the Newton method. The action can be
5688 /// overloaded at an elemental level by chaging the default
5689 /// behaviour of the function Element::get_all_vectors_and_matrices().
5690 /// column_or_row_index: Column [or row] index of given entry
5691 /// row_or_column_start: Index of first entry for given row [or column]
5692 /// value : Vector of nonzero entries
5693 /// residuals : Residual vector
5694 /// compressed_row_flag: Bool flag to indicate if storage format is
5695 /// compressed row [if false interpretation of
5696 /// arguments is as stated in square brackets].
5697 //=====================================================================
5699 Vector<int*>& column_or_row_index,
5700 Vector<int*>& row_or_column_start,
5701 Vector<double*>& value,
5702 Vector<unsigned>& nnz,
5703 Vector<double*>& residuals,
5704 bool compressed_row_flag)
5705 {
5706 // Total number of elements
5707 const unsigned long n_elements = mesh_pt()->nelement();
5708
5709 // Default range of elements for distributed problems
5710 unsigned long el_lo = 0;
5711 unsigned long el_hi = n_elements - 1;
5712
5713
5714#ifdef OOMPH_HAS_MPI
5715 // Otherwise just loop over a fraction of the elements
5716 // (This will either have been initialised in
5717 // Problem::set_default_first_and_last_element_for_assembly() or
5718 // will have been re-assigned during a previous assembly loop
5719 // Note that following the re-assignment only the entries
5720 // for the current processor are relevant.
5722 {
5723 el_lo = First_el_for_assembly[Communicator_pt->my_rank()];
5724 el_hi = Last_el_plus_one_for_assembly[Communicator_pt->my_rank()] - 1;
5725 }
5726#endif
5727
5728 // number of local eqns
5729 unsigned ndof = this->ndof();
5730
5731 // Find the number of vectors to be assembled
5732 const unsigned n_vector = residuals.size();
5733
5734 // Find the number of matrices to be assembled
5735 const unsigned n_matrix = column_or_row_index.size();
5736
5737 // Locally cache pointer to assembly handler
5739
5740#ifdef OOMPH_HAS_MPI
5741 bool doing_residuals = false;
5742 if (dynamic_cast<ParallelResidualsHandler*>(Assembly_handler_pt) != 0)
5743 {
5744 doing_residuals = true;
5745 }
5746#endif
5747
5748// Error check dimensions
5749#ifdef PARANOID
5750 if (row_or_column_start.size() != n_matrix)
5751 {
5752 std::ostringstream error_stream;
5753 error_stream << "Error: " << std::endl
5754 << "row_or_column_start.size() "
5755 << row_or_column_start.size() << " does not equal "
5756 << "column_or_row_index.size() "
5757 << column_or_row_index.size() << std::endl;
5758 throw OomphLibError(
5759 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5760 }
5761
5762 if (value.size() != n_matrix)
5763 {
5764 std::ostringstream error_stream;
5765 error_stream << "Error: " << std::endl
5766 << "value.size() " << value.size() << " does not equal "
5767 << "column_or_row_index.size() "
5768 << column_or_row_index.size() << std::endl
5769 << std::endl
5770 << std::endl;
5771 throw OomphLibError(
5772 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5773 }
5774#endif
5775
5776 // The idea behind this sparse assembly routine is to use Vectors of
5777 // Vectors for the entries in each complete matrix. And a second
5778 // Vector of Vectors stores the global row (or column) indeces. This
5779 // will not have the memory overheads associated with the methods using
5780 // lists or maps, but insertion will be more costly.
5781
5782 // Set up two vector of vectors to store the entries of each matrix,
5783 // indexed by either the column or row. The entries of the vector for
5784 // each matrix correspond to all the rows or columns of that matrix.
5785 Vector<Vector<Vector<unsigned>>> matrix_row_or_col_indices(n_matrix);
5786 Vector<Vector<Vector<double>>> matrix_values(n_matrix);
5787
5788 // Loop over the number of matrices being assembled and resize
5789 // each vector of vectors to the number of rows or columns of the matrix
5790 for (unsigned m = 0; m < n_matrix; m++)
5791 {
5792 matrix_row_or_col_indices[m].resize(ndof);
5793 matrix_values[m].resize(ndof);
5794 }
5795
5796 // Resize the residuals vectors
5797 for (unsigned v = 0; v < n_vector; v++)
5798 {
5799 residuals[v] = new double[ndof];
5800 for (unsigned i = 0; i < ndof; i++)
5801 {
5802 residuals[v][i] = 0;
5803 }
5804 }
5805
5806#ifdef OOMPH_HAS_MPI
5807
5808 // Storage for assembly time for elements
5809 double t_assemble_start = 0.0;
5810
5811 // Storage for assembly times
5812 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5813 {
5814 Elemental_assembly_time.resize(n_elements);
5815 }
5816
5817#endif
5818
5819
5820 //----------------Assemble and populate the vector storage scheme-------
5821 {
5822 // Allocate local storage for the element's contribution to the
5823 // residuals vectors and system matrices of the size of the maximum
5824 // number of dofs in any element
5825 // This means that the storage will only be allocated (and deleted) once
5826 Vector<Vector<double>> el_residuals(n_vector);
5827 Vector<DenseMatrix<double>> el_jacobian(n_matrix);
5828
5829 // Loop over the elements
5830 for (unsigned long e = el_lo; e <= el_hi; e++)
5831 {
5832#ifdef OOMPH_HAS_MPI
5833 // Time it?
5834 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5835 {
5836 t_assemble_start = TimingHelpers::timer();
5837 }
5838#endif
5839
5840 // Get the pointer to the element
5841 GeneralisedElement* elem_pt = mesh_pt()->element_pt(e);
5842
5843#ifdef OOMPH_HAS_MPI
5844 // Ignore halo elements
5845 if (!elem_pt->is_halo())
5846 {
5847#endif
5848
5849 // Find number of degrees of freedom in the element
5850 const unsigned nvar = assembly_handler_pt->ndof(elem_pt);
5851
5852 // Resize the storage for elemental jacobian and residuals
5853 for (unsigned v = 0; v < n_vector; v++)
5854 {
5855 el_residuals[v].resize(nvar);
5856 }
5857 for (unsigned m = 0; m < n_matrix; m++)
5858 {
5859 el_jacobian[m].resize(nvar);
5860 }
5861
5862 // Now get the residuals and jacobian for the element
5864 elem_pt, el_residuals, el_jacobian);
5865
5866 //---------------Insert the values into the vectors--------------
5867
5868 // Loop over the first index of local variables
5869 for (unsigned i = 0; i < nvar; i++)
5870 {
5871 // Get the local equation number
5872 unsigned eqn_number = assembly_handler_pt->eqn_number(elem_pt, i);
5873
5874 // Add the contribution to the residuals
5875 for (unsigned v = 0; v < n_vector; v++)
5876 {
5877 // Fill in each residuals vector
5878 residuals[v][eqn_number] += el_residuals[v][i];
5879 }
5880
5881 // Now loop over the other index
5882 for (unsigned j = 0; j < nvar; j++)
5883 {
5884 // Get the number of the unknown
5885 unsigned unknown = assembly_handler_pt->eqn_number(elem_pt, j);
5886
5887 // Loop over the matrices
5888 // If it's compressed row storage, then our vector of maps
5889 // is indexed by row (equation number)
5890 for (unsigned m = 0; m < n_matrix; m++)
5891 {
5892 // Get the value of the matrix at this point
5893 double value = el_jacobian[m](i, j);
5894 // Only bother to add to the vector if it's non-zero
5895 if (std::fabs(value) > Numerical_zero_for_sparse_assembly)
5896 {
5897 // If it's compressed row storage, then our vector of maps
5898 // is indexed by row (equation number)
5899 if (compressed_row_flag)
5900 {
5901 // Find the correct position and add the data into the
5902 // vectors
5903 const unsigned size =
5904 matrix_row_or_col_indices[m][eqn_number].size();
5905
5906 for (unsigned k = 0; k <= size; k++)
5907 {
5908 if (k == size)
5909 {
5910 matrix_row_or_col_indices[m][eqn_number].push_back(
5911 unknown);
5912 matrix_values[m][eqn_number].push_back(value);
5913 break;
5914 }
5915 else if (matrix_row_or_col_indices[m][eqn_number][k] ==
5916 unknown)
5917 {
5918 matrix_values[m][eqn_number][k] += value;
5919 break;
5920 }
5921 }
5922 }
5923 // Otherwise it's compressed column storage and our vector is
5924 // indexed by column (the unknown)
5925 else
5926 {
5927 // Add the data into the vectors in the correct position
5928 const unsigned size =
5929 matrix_row_or_col_indices[m][unknown].size();
5930 for (unsigned k = 0; k <= size; k++)
5931 {
5932 if (k == size)
5933 {
5934 matrix_row_or_col_indices[m][unknown].push_back(
5935 eqn_number);
5936 matrix_values[m][unknown].push_back(value);
5937 break;
5938 }
5939 else if (matrix_row_or_col_indices[m][unknown][k] ==
5940 eqn_number)
5941 {
5942 matrix_values[m][unknown][k] += value;
5943 break;
5944 }
5945 }
5946 }
5947 }
5948 } // End of loop over matrices
5949 }
5950 }
5951
5952#ifdef OOMPH_HAS_MPI
5953 } // endif halo element
5954#endif
5955
5956
5957#ifdef OOMPH_HAS_MPI
5958 // Time it?
5959 if ((!doing_residuals) && Must_recompute_load_balance_for_assembly)
5960 {
5962 TimingHelpers::timer() - t_assemble_start;
5963 }
5964#endif
5965
5966 } // End of loop over the elements
5967