42#include "triangle_mesh.h"
74 Explicit_time_stepper_pt(0),
76 Default_set_initial_condition_called(
false),
77 Use_globally_convergent_newton_method(
false),
78 Empty_actions_before_read_unstructured_meshes_has_been_called(
false),
79 Empty_actions_after_read_unstructured_meshes_has_been_called(
false),
80 Store_local_dof_pt_in_elements(
false),
81 Calculate_hessian_products_analytic(
false),
83 Doc_imbalance_in_parallel_assembly(
false),
84 Use_default_partition_in_load_balance(
false),
85 Must_recompute_load_balance_for_assembly(
true),
88 Relaxation_factor(1.0),
89 Newton_solver_tolerance(1.0e-8),
90 Max_newton_iterations(10),
91 Nnewton_iter_taken(0),
93 Time_adaptive_newton_crash_on_solve_fail(
false),
94 Jacobian_reuse_is_enabled(
false),
95 Jacobian_has_been_computed(
false),
96 Problem_is_nonlinear(
true),
97 Pause_at_end_of_sparse_assembly(
false),
98 Doc_time_in_distribute(
false),
99 Sparse_assembly_method(Perform_assembly_using_vectors_of_pairs),
100 Sparse_assemble_with_arrays_initial_allocation(400),
101 Sparse_assemble_with_arrays_allocation_increment(150),
102 Numerical_zero_for_sparse_assembly(0.0),
103 FD_step_used_in_get_hessian_vector_products(1.0e-8),
104 Mass_matrix_reuse_is_enabled(
false),
105 Mass_matrix_has_been_computed(
false),
106 Discontinuous_element_formulation(
false),
109 DTSF_max_increase(4.0),
110 DTSF_min_decrease(0.8),
111 Target_error_safety_factor(1.0),
112 Minimum_dt_but_still_proceed(-1.0),
113 Scale_arc_length(
true),
114 Desired_proportion_of_arc_length(0.5),
117 Continuation_direction(1.0),
118 Parameter_derivative(1.0),
119 Parameter_current(0.0),
120 Use_continuation_timestepper(
false),
121 Dof_derivative_offset(1),
122 Dof_current_offset(2),
124 Desired_newton_iterations_ds(5),
126 Bifurcation_detection(
false),
127 Bisect_to_find_bifurcation(
false),
128 First_jacobian_sign_change(
false),
129 Arc_length_step_taken(
false),
130 Use_finite_differences_for_continuation_derivatives(
false),
132 Dist_problem_matrix_distribution(Uniform_matrix_distribution),
133 Parallel_sparse_assemble_previous_allocation(0),
134 Problem_has_been_distributed(
false),
135 Bypass_increase_in_dof_check_during_pruning(
false),
136 Max_permitted_error_for_halo_check(1.0e-14),
138 Shut_up_in_newton_solve(
false),
139 Always_take_one_newton_step(
false),
140 Timestep_reduction_factor_after_nonconvergence(0.5),
141 Keep_temporal_error_below_tolerance(
true)
155#if defined(OOMPH_HAS_MUMPS) && \
156 defined(OOMPH_ENABLE_MUMPS_AS_DEFAULT_LINEAR_SOLVER)
217 for (
unsigned c = 0; c <
n_copies; c++)
266 for (
unsigned n = 0;
n <
n_var;
n++)
277 unsigned Nelement = 0;
313 const unsigned nrow = this->
ndof();
353 for (
unsigned p = 0;
p <
nproc;
p++)
377 error_stream <<
"Never get here. Dist_problem_matrix_distribution = "
451 for (
unsigned i = 0;
i <
n;
i++)
462 warn_message <<
"WARNING: All entries in specified partitioning vector \n"
463 <<
" are zero -- will ignore this and use METIS\n"
464 <<
" to perform the partitioning\n";
519 warn_message <<
"WARNING: You've tried to distribute a problem over\n"
520 <<
"only one processor: this would make METIS crash.\n"
521 <<
"Ignoring your request for distribution.\n";
523 "Problem::distribute()",
531 error_stream <<
"You have tried to distribute a problem\n"
532 <<
"but there are less elements than processors.\n"
533 <<
"Please re-run with more elements!\n"
534 <<
"Please also ensure that actions_before_distribute().\n"
535 <<
"and actions_after_distribute() are correctly set up.\n"
590 error_stream <<
"You have tried to distribute a problem\n"
591 <<
"but at least one of your meshes is no longer\n"
592 <<
"uniformly refined. In order to preserve the Tree\n"
593 <<
"and TreeForest structure, Problem::distribute() can\n"
594 <<
"only be called while meshes are uniformly refined.\n"
605 error_stream <<
"You have tried to distribute a problem\n"
606 <<
"and there is some global data.\n"
607 <<
"This is not likely to work...\n"
653 oomph_info <<
"INFO: using METIS to partition elements" << std::endl;
659 oomph_info <<
"INFO: using pre-set partition of elements"
672 oomph_info <<
"Time for partitioning of global mesh: "
696 oomph_info <<
"Time for actions before distribute: "
765 for (
unsigned e = 0;
e <
n_ele;
e++)
815 for (
unsigned e = 0;
e <
n_ele;
e++)
832 <<
") is not the sameas the number of\nadded elements ("
833 <<
counter <<
") to the Base_mesh_element_pt data "
834 <<
"structure!!!\n\n";
836 "Problem::distribute()",
851 oomph_info <<
"INFO: We're over-ruling the \"keep as halo element\"\n"
852 <<
" status because the preset partitioning\n"
853 <<
" didn't place ANY elements on this processor,\n"
854 <<
" probably because of a restart on a larger \n"
855 <<
" number of processors\n";
878 <<
"--------------------" << std::endl;
883 submesh_element_domain[
i_mesh],
895 for (
unsigned e = 0;
e <
n_del;
e++)
932 oomph_info <<
"Time for re-assigning eqn numbers (in distribute): "
942 <<
"Number of dofs in distribute() has changed "
944 <<
"Check that you've implemented any necessary "
945 "actions_before/after\n"
946 <<
"distribute functions, e.g. to pin redundant pressure dofs"
993 << doc_info.
number() <<
".dat";
1024 for (
unsigned e = 0;
e <
nelem;
e++)
1033 for (
unsigned e = 0;
e <
nelem;
e++)
1056 for (
int d = 0; d <
n_proc; d++)
1078 error_stream <<
"No process has more than 1 element, and\n"
1079 <<
"at least one process has no elements!\n"
1080 <<
"Suggest rerunning with more refinement.\n"
1089 for (
unsigned e = 0;
e <
nelem;
e++)
1102 oomph_info <<
"INFO: Switched element domain at position " <<
e
1105 <<
" to process " << d << std::endl
1106 <<
"which was given no elements by METIS partition"
1119 for (
unsigned e = 0;
e <
nelem;
e++)
1125 for (
unsigned e = 0;
e <
nelem;
e++)
1136 <<
" elements from this partition" << std::endl
1159 <<
"WARNING: Problem::prune_halo_elements_and_nodes() was called on a "
1160 <<
"non-distributed Problem!" << std::endl;
1161 oomph_info <<
"Ignoring your request..." << std::endl;
1169 <<
"WARNING: You've tried to prune halo layers on a problem\n"
1170 <<
"with only one processor: this is unnecessary.\n"
1171 <<
"Ignoring your request." << std::endl
1193 oomph_info <<
"Time for actions_before_distribute() in "
1194 <<
"Problem::prune_halo_elements_and_nodes(): "
1201 std::map<GeneralisedElement*, unsigned>
1204 for (
unsigned e = 0;
e <
nel;
e++)
1229 ref_el_pt->tree_pt()->stick_all_tree_nodes_into_vector(tree_pt);
1230 unsigned ntree = tree_pt.
size();
1231 for (
unsigned t = 0;
t < ntree;
t++)
1244 oomph_info <<
"Time for setup old root elements in "
1245 <<
"Problem::prune_halo_elements_and_nodes(): "
1280 oomph_info <<
"Total time for all mesh-level prunes in "
1281 <<
"Problem::prune_halo_elements_and_nodes(): "
1373 ref_el_pt->tree_pt()->root_pt()->object_pt();
1410 oomph_info <<
"Time for setup of new base elements in "
1411 <<
"Problem::prune_halo_elements_and_nodes(): "
1432 <<
"Problem::prune_halo_elements_and_nodes(): "
1457 <<
"Number of new root elements spawned from old root " <<
e
1517 oomph_info <<
"Time for finishing off base mesh info "
1518 <<
"Problem::prune_halo_elements_and_nodes(): "
1531 oomph_info <<
"Time for actions_after_distribute() "
1532 <<
"Problem::prune_halo_elements_and_nodes(): "
1549 oomph_info <<
"Time for assign_eqn_numbers() "
1550 <<
"Problem::prune_halo_elements_and_nodes(): "
1563 <<
"Number of dofs in prune_halo_elements_and_nodes() has "
1566 <<
"Check that you've implemented any necessary "
1567 "actions_before/after"
1568 <<
"\nadapt/distribute functions, e.g. to pin redundant pressure"
1595 std::string error_message =
"Problem::build_global_mesh() called,\n";
1596 error_message +=
" but a global mesh has already been built:\n";
1597 error_message +=
"Problem::Mesh_pt is not zero!\n";
1605 std::string error_message =
"Problem::build_global_mesh() called,\n";
1606 error_message +=
" but there are no submeshes:\n";
1607 error_message +=
"Problem::Sub_mesh_pt has no entries\n";
1654 oomph_info <<
"Created Time with " << ndt <<
" timesteps" << std::endl;
1663 oomph_info <<
"Resized Time to include " << ndt <<
" timesteps"
1669 oomph_info <<
"Time object already has storage for " << ndt
1670 <<
" timesteps" << std::endl;
1693 oomph_info <<
"Created Time with storage for no previous timestep"
1698 oomph_info <<
"Time object already exists " << std::endl;
1759 oomph_info <<
"Problem is not distributed. Parallel assembly of "
1760 <<
"Jacobian uses default partitioning: " << std::endl;
1765 oomph_info <<
"Proc " <<
p <<
" assembles from element "
1771 oomph_info <<
"Proc " <<
p <<
" assembles no elements\n";
1799 oomph_info <<
"Not re-computing distribution of elemental assembly\n"
1800 <<
"because there are fewer elements than processors\n";
1824 for (
unsigned e = 0;
e <
nel;
e++)
1857 <<
"Re-assigning distribution of element assembly over processors:"
1928 error_stream <<
"Error: First/last element of proc " <<
p <<
": "
1936 error_stream <<
"Error: First element of proc " <<
p + 1 <<
": "
2018 oomph_info <<
"Processor " << 0 <<
" assembles Jacobians"
2037 oomph_info <<
"Processor " <<
p <<
" assembles Jacobians"
2086 const bool& assign_local_eqn_numbers)
2093 error_stream <<
"Global mesh does not exist, so equation numbers cannot "
2098 error_stream <<
"There aren't even any sub-meshes in the Problem.\n"
2099 <<
"You can set the global mesh directly by using\n"
2100 <<
"Problem::mesh_pt() = my_mesh_pt;\n"
2101 <<
"OR you can use Problem::add_sub_mesh(mesh_pt); "
2102 <<
"to add a sub mesh.\n";
2108 error_stream <<
"You need to call Problem::build_global_mesh() to create "
2110 <<
"from the sub-meshes.\n\n";
2161 for (
unsigned e = 0;
e <
nel;
e++)
2190 <<
"Time for complete setup of dependencies in assign_eqn_numbers: "
2235 n_dof = spine_mesh_pt->assign_global_spine_eqn_numbers(
Dof_pt);
2247 n_dof = spine_mesh_pt->assign_global_spine_eqn_numbers(
Dof_pt);
2256 <<
"Time for assign_global_eqn_numbers in assign_eqn_numbers: "
2289 oomph_info <<
"Time for Problem::synchronise_eqn_numbers in "
2333 std::stringstream
tmp;
2334 tmp <<
"Time for calls to Problem::remove_duplicate_data in "
2341 tmp <<
" removed some/any data.\n";
2363 std::stringstream
tmp;
2365 <<
"Time for MPI_Allreduce after Problem::remove_duplicate_data in "
2366 <<
"Problem::assign_eqn_numbers: " <<
t_end -
t_start << std::endl;
2388 <<
"Problem::remove_null_pointers_from_external_halo_node_"
2420 if (assign_local_eqn_numbers)
2439 oomph_info <<
"Total time for all Mesh::assign_local_eqn_numbers in "
2463 <<
"Global mesh does not exist, so equation numbers cannot be found.\n";
2467 error_stream <<
"There aren't even any sub-meshes in the Problem.\n"
2468 <<
"You can set the global mesh directly by using\n"
2469 <<
"Problem::mesh_pt() = my_mesh_pt;\n"
2470 <<
"OR you can use Problem::add_sub_mesh(mesh_pt); "
2471 <<
"to add a sub mesh.\n";
2477 error_stream <<
"You need to call Problem::build_global_mesh() to create "
2479 <<
"from the sub-meshes.\n\n";
2487 <<
"Although this program will describe the degrees of freedom in the \n"
2488 <<
"problem, it will do so using the typedef for the elements. This is \n"
2489 <<
"not neccesarily human readable, but there is a solution.\n"
2490 <<
"Pipe your program's output through c++filt, with the argument -t.\n"
2491 <<
"e.g. \"./two_d_multi_poisson | c++filt -t > ReadableOutput.txt\".\n "
2492 <<
"(Disregarding the quotes)\n\n\n";
2494 out <<
"Classifying Global Equation Numbers" << std::endl;
2495 out << std::string(80,
'-') << std::endl;
2513 std::string
in(
" in Problem's Only Mesh.");
2523 std::string
in(
" in Problem's Only SpineMesh.");
2524 spine_mesh_pt->describe_spine_dofs(
out,
in);
2539 spine_mesh_pt->describe_spine_dofs(
out,
in);
2545 out << std::string(80,
'\\') << std::endl;
2546 out << std::string(80,
'\\') << std::endl;
2547 out << std::string(80,
'\\') << std::endl;
2548 out <<
"Classifying global eqn numbers in terms of elements." << std::endl;
2549 out << std::string(80,
'-') << std::endl;
2550 out <<
"Eqns | Source" << std::endl;
2551 out << std::string(80,
'-') << std::endl;
2555 std::string
in(
" in Problem's Only Mesh.");
2584 for (
unsigned long l = 0;
l <
n_dof;
l++)
2596 throw OomphLibError(
"Not designed for distributed problems",
2614 if (eqn_number >= 0)
2631 if (eqn_number >= 0)
2647 if (eqn_number >= 0)
2679 actually_removed_some_data =
false;
2770 ->ncont_interpolated_values();
2939 unsigned nb = (*boundaries_pt).
size();
2941 for (std::set<unsigned>::iterator
it =
2942 (*boundaries_pt).begin();
2943 it != (*boundaries_pt).end();
2948 for (
unsigned i = 0;
i <
nb;
i++)
2977 ->ncont_interpolated_values();
3004 <<
"Number of master nodes for node to be replaced, "
3007 <<
" for i_cont=" <<
i_cont << std::endl;
3010 <<
"Nodal coordinates of replacement node:";
3012 for (
unsigned i = 0;
i < ndim;
i++)
3019 <<
" master nodes are: \n";
3024 ->master_node_pt(
k);
3026 for (
unsigned i = 0;
i < ndim;
i++)
3036 <<
"Nodal coordinates of node to be replaced:";
3038 for (
unsigned i = 0;
i < ndim;
i++)
3046 <<
" master nodes are: \n";
3053 for (
unsigned i = 0;
i < ndim;
i++)
3088 ->ncont_interpolated_values();
3166 for (std::set<unsigned>::iterator
it =
3167 (*boundaries_pt).begin();
3168 it != (*boundaries_pt).end();
3195 <<
"About to re-set master for i_cont= " <<
i_cont
3196 <<
" for external node (with proc " <<
iproc
3199 for (
unsigned jj = 0;
jj <
n;
jj++)
3204 <<
" which is not hanging --> About to die!"
3205 <<
"Outputting offending element into oomph-info "
3342 for (
unsigned j = 0;
j <
nnod;
j++)
3479 for (
unsigned j = 0;
j <
nnod;
j++)
3509 const unsigned long n_dof = this->
ndof();
3511 if (n_dof !=
dofs.nrow())
3514 error_stream <<
"Number of degrees of freedom in vector argument "
3515 <<
dofs.nrow() <<
"\n"
3516 <<
"does not equal number of degrees of freedom in problem "
3522 for (
unsigned long l = 0;
l <
n_dof;
l++)
3534 throw OomphLibError(
"Not designed for distributed problems",
3549 if (eqn_number >= 0)
3566 if (eqn_number >= 0)
3582 if (eqn_number >= 0)
3598 throw OomphLibError(
"Not implemented for distributed problems!",
3616 if (eqn_number >= 0)
3632 if (eqn_number >= 0)
3648 if (eqn_number >= 0)
3663 const unsigned long n_dof = this->
ndof();
3664 for (
unsigned long l = 0;
l <
n_dof;
l++)
3685 error_stream <<
"The function get_inverse_mass_matrix_times_residuals() "
3687 <<
"used with the default assembly handler\n\n";
3734 oomph_info <<
"Not recomputing Mass Matrix " << std::endl;
3753 oomph_info <<
"Enabling resolve in explicit timestep" << std::endl;
3832 error_stream <<
"The distribution of the residuals vector does not "
3833 "have the correct\n"
3834 <<
"number of global rows\n";
3912 dist_pt, column_index, row_start, value, nnz,
res);
3947 if (
residuals.distribution_pt()->distributed())
3951 <<
"If the DoubleVector residuals is setup then it must not "
3952 <<
"be distributed.";
3960 <<
"If the DoubleVector residuals is setup then it must have"
3961 <<
" the correct number of rows";
3966 *
residuals.distribution_pt()->communicator_pt()))
3970 <<
"If the DoubleVector residuals is setup then it must have"
3971 <<
" the same communicator as the problem.";
4076 error_stream <<
"The distribution of the residuals must "
4077 <<
"be the same as the distribution of the jacobian.";
4085 <<
"The distribution of the jacobian and residuals does not"
4086 <<
"have the correct number of global rows.";
4094 error_stream <<
"The distribution of the jacobian and residuals must "
4095 <<
"both be setup or both not setup";
4132 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4142 dist_pt, column_index, row_start, value, nnz,
res);
4145 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4157 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4202 if (
residuals.distribution_pt()->distributed())
4206 <<
"If the DoubleVector residuals is setup then it must not "
4207 <<
"be distributed.";
4215 <<
"If the DoubleVector residuals is setup then it must have"
4216 <<
" the correct number of rows";
4221 *
residuals.distribution_pt()->communicator_pt()))
4225 <<
"If the DoubleVector residuals is setup then it must have"
4226 <<
" the same communicator as the problem.";
4269 value[0], row_index[0], column_start[0], nnz[0],
n_dof,
n_dof);
4277 error_stream <<
"Cannot assemble a CCDoubleMatrix Jacobian on more "
4278 <<
"than one processor.";
4349 for (
unsigned i = 0;
i <
n_dim;
i++)
4417 for (
unsigned i = 0;
i <
n_dim;
i++)
4549 <<
"Error: Incorrect value for Problem::Sparse_assembly_method"
4551 <<
"It should be one of the enumeration Problem::Assembly_method"
4587 unsigned long el_lo = 0;
4630 <<
"row_or_column_start.size() "
4632 <<
"column_or_row_index.size() "
4642 <<
"Error in Problem::sparse_assemble_row_or_column_compressed "
4644 <<
"value.size() " << value.
size() <<
" does not equal "
4685 for (
unsigned i = 0;
i <
ndof;
i++)
4755 for (
unsigned i = 0;
i <
nvar;
i++)
4768 for (
unsigned j = 0;
j <
nvar;
j++)
4863 value[
m] =
new double[nnz[
m]];
4879 for (std::map<unsigned, double>::iterator
it =
4899 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
4900 pause(
"Check memory usage now.");
4933 unsigned long el_lo = 0;
4976 <<
"row_or_column_start.size() "
4978 <<
"column_or_row_index.size() "
4988 <<
"Error in Problem::sparse_assemble_row_or_column_compressed "
4990 <<
"value.size() " << value.
size() <<
" does not equal "
5029 for (
unsigned i = 0;
i <
ndof;
i++)
5060 std::list<std::pair<unsigned, double>>*
list_pt;
5102 for (
unsigned i = 0;
i <
nvar;
i++)
5115 for (
unsigned j = 0;
j <
nvar;
j++)
5138 std::make_pair(
unknown, value));
5150 std::make_pair(eqn_number, value));
5219 value[
m] =
new double[nnz[
m]];
5241 std::list<std::pair<unsigned, double>>::iterator
it =
5314 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
5315 pause(
"Check memory usage now.");
5348 unsigned long el_lo = 0;
5391 <<
"row_or_column_start.size() "
5393 <<
"column_or_row_index.size() "
5403 <<
"value.size() " << value.
size() <<
" does not equal "
5404 <<
"column_or_row_index.size() "
5434 for (
unsigned i = 0;
i <
ndof;
i++)
5502 for (
unsigned i = 0;
i <
nvar;
i++)
5515 for (
unsigned j = 0;
j <
nvar;
j++)
5537 for (
unsigned k = 0;
k <= size;
k++)
5542 std::make_pair(
unknown, value));
5558 for (
unsigned k = 0;
k <= size;
k++)
5563 std::make_pair(eqn_number, value));
5630 for (
unsigned long i = 0;
i <
ndof;
i++)
5667 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
5668 pause(
"Check memory usage now.");
5701 unsigned long el_lo = 0;
5745 <<
"row_or_column_start.size() "
5747 <<
"column_or_row_index.size() "
5757 <<
"value.size() " << value.
size() <<
" does not equal "
5758 <<
"column_or_row_index.size() "
5791 for (
unsigned i = 0;
i <
ndof;
i++)
5860 for (
unsigned i = 0;
i <
nvar;
i++)
5873 for (
unsigned j = 0;
j <
nvar;
j++)
5894 const unsigned size =
5897 for (
unsigned k = 0;
k <= size;
k++)
5919 const unsigned size =
5921 for (
unsigned k = 0;
k <= size;
k++)
5993 for (
unsigned long i = 0;
i <
ndof;
i++)
6030 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
6031 pause(
"Check memory usage now.");
6064 unsigned long el_lo = 0;
6108 <<
"row_or_column_start.size() "
6110 <<
"column_or_row_index.size() "
6120 <<
"value.size() " << value.
size() <<
" does not equal "
6121 <<
"column_or_row_index.size() "
6154 for (
unsigned i = 0;
i <
ndof;
i++)
6238 for (
unsigned i = 0;
i <
nvar;
i++)
6251 for (
unsigned j = 0;
j <
nvar;
j++)
6267 const unsigned size =
ncoef[
m][eqn_number];
6274 [
m][eqn_number] != 0)
6300 for (
unsigned k = 0;
k <= size;
k++)
6306 [
m][eqn_number] ==
ncoef[
m][eqn_number])
6313 for (
unsigned c = 0; c <
ncoef[
m][eqn_number]; c++)
6328 unsigned entry =
ncoef[
m][eqn_number];
6348 for (
unsigned k = 0;
k <= size;
k++)
6443 for (
unsigned long i = 0;
i <
ndof;
i++)
6488 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
6489 pause(
"Check memory usage now.");
6501 const unsigned&
el_lo,
6502 const unsigned&
el_hi,
6523 for (
unsigned i = 0;
i <
nvar;
i++)
6526 unsigned global_eqn_number =
6574 <<
"This is usually a sign that the problem distribution \n"
6575 <<
"or the load balancing have gone wrong.";
6577 "Problem::parallel_sparse_assemble()",
6584 unsigned long el_lo = 0;
6620 <<
"row_or_column_start.size() " << row_start.
size()
6621 <<
" does not equal "
6632 <<
"value.size() " << values.
size() <<
" does not equal "
6769 for (
unsigned i = 0;
i <
nvar;
i++)
6772 unsigned global_eqn_number =
6779 int eqn_number =
right / 2;
6780 while (
my_eqns[eqn_number] != global_eqn_number)
6798 for (
unsigned j = 0;
j <
nvar;
j++)
6824 <<
"Internal Error: " << std::endl
6825 <<
"Could not find global equation number "
6826 << global_eqn_number
6827 <<
" in my_eqns vector of equation numbers but\n"
6828 <<
"at least one entry in the residual vector is nonzero.";
6838 if (
my_eqns[eqn_number] > global_eqn_number)
6857 for (
unsigned j = 0;
j <
nvar;
j++)
6873 const unsigned size =
ncoef[
m][eqn_number];
6880 [
m][eqn_number] != 0)
6905 for (
unsigned k = 0;
k <= size;
k++)
6911 [
m][eqn_number] ==
ncoef[
m][eqn_number])
6918 for (
unsigned c = 0; c <
ncoef[
m][eqn_number]; c++)
6933 unsigned entry =
ncoef[
m][eqn_number];
6996 oomph_info <<
"\nCPU for residual computation (loc/max/min/imbal): ";
7000 oomph_info <<
"\nCPU for Jacobian computation (loc/max/min/imbal): ";
7027 for (
unsigned c = 0; c <
ncoef[
m][
e]; c++)
7086 for (
unsigned p = 0;
p <
nproc;
p++)
7104 for (
unsigned p = 0;
p <
nproc;
p++)
7147 for (
unsigned p = 0;
p <
nproc;
p++)
7165 for (
unsigned p = 0;
p <
nproc;
p++)
7185 for (
unsigned p = 0;
p <
nproc;
p++)
7231 for (
unsigned p = 0;
p <
nproc;
p++)
7271 for (
unsigned p = 0;
p <
nproc;
p++)
7470 row_start[
m][0] = 0;
7475 for (
unsigned p = 0;
p <
nproc;
p++)
7491 row_start[
m][
i] = 0;
7495 for (
unsigned p = 0;
p <
nproc;
p++)
7550 for (
unsigned p = 0;
p <
nproc;
p++)
7575 for (
unsigned c = 0; c <
n_chunk; c++)
7642 for (
unsigned p = 0;
p <
nproc;
p++)
7657 for (
unsigned c = 0; c <
n_chunk; c++)
7664 values[
m] =
new double[nnz[
m]];
7669 for (
unsigned c = 0; c <
n_chunk; c++)
7672 for (
unsigned i = 0;
i <
nc;
i++)
7685 unsigned g = row_start[
m][0];
7686 row_start[
m][0] = 0;
7689 unsigned h = g + row_start[
m][
i];
7690 row_start[
m][
i] = g;
7704 for (
unsigned p = 0;
p <
nproc;
p++)
7719 for (
unsigned p = 0;
p <
nproc;
p++)
7730 for (
unsigned p = 0;
p <
nproc;
p++)
7746 for (
unsigned p = 0;
p <
nproc;
p++)
7797 oomph_info <<
"CPU for residual distribut. (loc/max/min/imbal): ";
7801 oomph_info <<
"CPU for Jacobian distribut. (loc/max/min/imbal): ";
7821 OomphLibWarning(
"This is unlikely to work with a distributed problem",
7822 " Problem::get_fd_jacobian()",
7837 const double FD_step = 1.0e-8;
7933 const double FD_step = 1.0e-8;
7993 for (
unsigned i = 0;
i <
n_vec;
i++)
8003 for (
unsigned i = 0;
i <
n_vec;
i++)
8038 for (
unsigned l = 0;
l <
n_var;
l++)
8041 const unsigned long eqn_number =
8045 for (
unsigned i = 0;
i <
n_vec;
i++)
8047 C_local(
i,
l) = C[
i].global_value(eqn_number);
8056 for (
unsigned l = 0;
l <
n_var;
l++)
8058 const unsigned long eqn_number =
8061 for (
unsigned i = 0;
i <
n_vec;
i++)
8097 for (
unsigned i = 0;
i <
n_vec;
i++)
8115 for (
unsigned i = 0;
i <
n_vec;
i++)
8131 for (
unsigned i = 0;
i <
n_vec;
i++)
8140 for (
unsigned i = 0;
i <
n_vec;
i++)
8154 for (
unsigned long e = 0;
e < n_element;
e++)
8173 for (
unsigned n = 0;
n <
n_var;
n++)
8180 for (
unsigned i = 0;
i <
n_vec;
i++)
8183 for (
unsigned n = 0;
n <
n_var;
n++)
8188 C_mult[
i] * C[
i].global_value(eqn_number);
8199 for (
unsigned n = 0;
n <
n_var;
n++)
8207 for (
unsigned n = 0;
n <
n_var;
n++)
8211 for (
unsigned m = 0;
m <
n_var;
m++)
8233 for (
unsigned i = 0;
i <
n_vec;
i++)
8235 product[
i].sum_all_halo_and_haloed_values();
8246 Vector<std::complex<double>>& alpha,
8313 Vector<std::complex<double>>& eigenvalue,
8378 Vector<std::complex<double>>& eigenvalue,
8444 const double&
shift)
8460 <<
"The distributions of the jacobian and mass matrix are\n"
8461 <<
"not the same and they must be.\n";
8470 <<
"mass_matrix has a distribution, but the number of rows is not "
8471 <<
"equal to the number of degrees of freedom in the problem.";
8480 <<
"main_matrix has a distribution, but the number of rows is not "
8481 <<
"equal to the number of degrees of freedom in the problem.";
8491 error_stream <<
"The distribution of the jacobian and mass matrix must "
8492 <<
"both be setup or both not setup";
8649 (*Saved_dof_pt)[
i] = *(this->
Dof_pt[
i]);
8663 for (
unsigned long n = 0;
n <
n_dof;
n++)
8665 (*Saved_dof_pt)[
n] =
dof(
n);
8679 "There are no stored values, use store_current_dof_values()\n",
8694 throw OomphLibError(
"The number of stored values is not equal to the "
8695 "current number of dofs\n",
8715 throw OomphLibError(
"The number of stored values is not equal to the "
8716 "current number of dofs\n",
8722 for (
unsigned long n = 0;
n <
n_dof;
n++)
8724 dof(
n) = (*Saved_dof_pt)[
n];
8742 std::ostringstream error_message;
8743 error_message <<
"Eigenvector has size " <<
eigenvector.nrow()
8744 <<
", not equal to the number of dofs in the problem,"
8745 <<
n_dof << std::endl;
8780 std::ostringstream error_message;
8781 error_message <<
"Eigenvector has size " <<
eigenvector.nrow()
8782 <<
", not equal to the number of dofs in the problem,"
8783 <<
n_dof << std::endl;
8852 error_stream <<
"Globally convergent Newton method has not been "
8853 <<
"implemented in parallel yet!" << std::endl;
8884 oomph_info << std::endl << std::endl << std::endl;
8885 oomph_info <<
"This is a bit bizarre: The problem has no dofs."
8888 <<
"I'll just return from the Newton solver without doing anything."
8897 oomph_info <<
"I hope this is what you intended me to do..."
8900 <<
"Note: All actions_...() functions were called"
8902 oomph_info << std::endl <<
" before returning." << std::endl;
8903 oomph_info << std::endl << std::endl << std::endl;
8939 double maxres =
dx.max();
8953 oomph_info <<
"\nInitial Maximum residuals " << maxres << std::endl;
8968 <<
"Linear problem -- convergence in one iteration assumed."
8986 oomph_info <<
"Not recomputing Jacobian! " << std::endl;
9006 oomph_info <<
"Enabling resolve" << std::endl;
9030 double*
dx_pt =
dx.values_pt();
9077 double maxres = 0.0;
9098 << maxres << std::endl
9119 <<
") has been exceeded in Newton solver." << std::endl;
9123 oomph_info <<
"Reached max. number of iterations ("
9159 oomph_info <<
"Time outside linear solver : "
9162 <<
" %" << std::endl;
9169 oomph_info <<
"Time outside linear solver : "
9170 <<
"[too fast]" << std::endl;
9194 for (
unsigned i = 0;
i <
n_dof;
i++)
9201 for (
unsigned i = 0;
i <
n_dof;
i++)
9207 for (
unsigned i = 0;
i <
n_dof;
i++)
9214 warn_message <<
"WARNING: Non-negative slope, probably due to a "
9215 <<
" roundoff \nproblem in the linesearch: slope=" <<
slope
9218 "Problem::globally_convergent_line_search()",
9222 for (
unsigned i = 0;
i <
n_dof;
i++)
9229 double lambda = 1.0;
9232 for (
unsigned i = 0;
i <
n_dof;
i++)
9241 for (
unsigned i = 0;
i <
n_dof;
i++)
9252 warn_message <<
"WARNING: Line search converged on x only!\n";
9254 "Problem::globally_convergent_line_search()",
9261 oomph_info <<
"Returning from linesearch with lambda=" << lambda
9357 <<
"USER-DEFINED ERROR IN NEWTON SOLVER " << std::endl;
9359 if (
error.linear_solver_error())
9361 oomph_info <<
"ERROR IN THE LINEAR SOLVER" << std::endl;
9367 <<
") REACHED WITHOUT CONVERGENCE " << std::endl;
9379 error_stream <<
"Error occured in Newton solver. " << std::endl;
9484 double maxres = y.
max();
9525 oomph_info <<
"Initial Maximum residuals " << maxres << std::endl;
9555 ->solve_for_two_rhs(
this, y,
input_z, z);
9662 double maxres = y.
max();
9902 if (n_sub_mesh == 0)
9906 if (spine_mesh_pt->does_pointer_correspond_to_spine_data(
parameter_pt))
9921 if (spine_mesh_pt->does_pointer_correspond_to_spine_data(
10056#ifdef OOMPH_HAS_MPI
10096 <<
"Warning: This function is called after spatially adapting the\n"
10097 <<
"eigenfunction associated with a pitchfork bifurcation and should\n"
10098 <<
"ensure that the exact (anti-)symmetries of problem are enforced\n"
10099 <<
"within that eigenfunction. It is problem specific and must be\n"
10100 <<
"filled in by the user if required.\n"
10101 <<
"A sign of problems is if the slack paramter gets too large and\n"
10102 <<
"if the solution at the Pitchfork is not symmetric.\n";
10105 "Problem::symmetrise_eigenfunction_for_adaptive_pitchfork_tracking()",
10280 const double& omega,
10335 std::ostringstream error_message;
10337 <<
"The parameter addressed by " <<
parameter_pt <<
" with the value "
10339 <<
"\n is supposed to be used for arc-length contiunation,\n"
10340 <<
" but it is stored in a Data object used by the problem.\n\n"
10341 <<
"This is bad for two reasons:\n"
10342 <<
"1. If it's a variable in the problem, it must already have an\n"
10343 "associated equation, so it can't be used for continuation;\n"
10344 <<
"2. The problem data will be reorganised in memory during "
10346 <<
" which means that the pointer will become invalid.\n\n"
10347 <<
"If you are sure that this is what you want to do you must:\n"
10348 <<
"A. Ensure that the value is pinned (don't worry we'll shout again "
10350 <<
"B. Use the alternative interface\n"
10351 <<
" Problem::arc_length_step_solve(Data*,unsigned,...)\n"
10352 <<
" which uses a pointer to the data object and not the raw double "
10366 for (
unsigned i = 0;
i < n_time_steppers;
i++)
10370 continuation_time_stepper_added =
true;
10378 oomph_info <<
"Adding the continuation time stepper\n";
10394 <<
" equation numbers allocated for continuation\n";
10425 <<
" in the data object to be used for continuation\n"
10426 <<
"is not pinned, which means that it is already a\n"
10427 <<
"variable in the problem "
10428 <<
"and cannot be used for continuation.\n\n"
10429 <<
"Please correct your formulation by either:\n"
10430 <<
"A. Pinning the value"
10432 <<
"B. Using a different parameter for continuation"
10445 for (
unsigned i = 0;
i < n_time_steppers;
i++)
10449 continuation_time_stepper_added =
true;
10457 oomph_info <<
"Adding the continuation time stepper\n";
10474 <<
" equation numbers allocated for continuation\n";
10508 if (n_sub_mesh == 0)
10512 spine_mesh_pt->set_consistent_pinned_spine_values_for_continuation(
10526 spine_mesh_pt->set_consistent_pinned_spine_values_for_continuation(
10603#ifdef OOMPH_HAS_MPI
10629#ifdef OOMPH_HAS_MPI
10630 <<
", in total (over all processors).\n";
10639 oomph_info <<
"\n \n Solution is fully converged in "
10640 <<
"Problem::newton_solver(). \n \n ";
10686 unsigned count = 0;
10704 std::ostringstream error_message;
10705 error_message <<
"DESIRED ARC-LENGTH STEP " <<
Ds_current
10706 <<
" HAS FALLEN BELOW MINIMUM TOLERANCE, " <<
Minimum_ds
10740 if (
error.linear_solver_error())
10744 <<
"USER-DEFINED ERROR IN NEWTON SOLVER " << std::endl;
10745 oomph_info <<
"ERROR IN THE LINEAR SOLVER" << std::endl;
10753 oomph_info <<
"STEP REJECTED DUE TO NEWTON SOLVER --- TRYING AGAIN"
10763 <<
"STEP REJECTED DUE TO INVERTED ELEMENTS --- TRYING AGAIN"
10788 std::string error_message =
10789 "The sign of the jacobian is zero after a linear solve\n";
10790 error_message +=
"Either the matrix is singular (unlikely),\n";
10791 error_message +=
"or the linear solver cannot compute the "
10792 "determinant of the matrix;\n";
10793 error_message +=
"e.g. an iterative linear solver.\n";
10795 "If the latter, bifurcation detection must be via an eigensolver\n";
10797 "Problem::arc_length_step_solve",
10832 <<
"-----------------------------------------------------------";
10834 <<
"SIGN CHANGE IN DETERMINANT OF JACOBIAN: " << std::endl;
10835 message <<
"BIFURCATION OR TURNING POINT DETECTED BETWEEN "
10841 <<
"-----------------------------------------------------------"
10849 std::ios_base::app);
10966 throw OomphLibError(
"Explicit time stepper pointer is null in problem.",
11052 <<
"USER-DEFINED ERROR IN NEWTON SOLVER " << std::endl;
11054 if (
error.linear_solver_error())
11056 oomph_info <<
"ERROR IN THE LINEAR SOLVER" << std::endl;
11062 <<
") REACHED WITHOUT CONVERGENCE " << std::endl;
11073 error_stream <<
"Error occured in unsteady Newton solver. " << std::endl;
11100 const double& epsilon)
11120 const double& epsilon,
11142 unsigned adaptive_flag = 0;
11209 if (
error.linear_solver_error() ||
11212 std::string error_message =
"USER-DEFINED ERROR IN NEWTON SOLVER\n";
11213 error_message +=
"ERROR IN THE LINEAR SOLVER\n";
11222 oomph_info <<
"TIMESTEP REJECTED DUE TO THE NEWTON SOLVER"
11233 oomph_info <<
"TIMESTEP REJECTED DUE TO INVERTED ELEMENTS" << std::endl;
11278 <<
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n"
11279 <<
"Estimated timestepping error is " <<
error <<
"\n"
11284 if (
error > epsilon)
11287 <<
" exceeds tolerance " << epsilon <<
"\n";
11295 oomph_info <<
" ...but we're not rejecting the timestep\n";
11298 <<
"Note: This behaviour can be adjusted by changing the\n"
11299 <<
"protected boolean\n"
11300 <<
" Problem::Keep_temporal_error_below_tolerance\n\n"
11301 <<
"Also, if you are noticing that many of your timesteps result\n"
11302 <<
"in error > tolerance, try reducing the target error with\n"
11303 <<
"respect to the error tolerance by reducing the value of\n"
11304 <<
"Target_error_safety_factor from its default value of 1.0\n"
11305 <<
"using the access function\n"
11306 <<
" target_error_safety_factor() = 0.5 (e.g.)\n"
11307 <<
"The default strategy (Target_error_safety_factor=1.0) tries\n"
11308 <<
"to suggest a timestep which will produce an error equal to\n"
11309 <<
"the error tolerance `epsilon` which risks error > tolerance\n"
11310 <<
"quite often. Setting the safety factor to too small a value\n"
11311 <<
"will make the timesteps unnecessarily small; too large will\n"
11312 <<
"not address the issue -- neither is optimal and a problem\n"
11313 <<
"dependent compromise is needed.\n"
11314 <<
"for more info see:\n"
11315 <<
" Mayr et al. (2018), p5,9, DOI:10.1016/j.finel.2017.12.002\n"
11316 <<
" Harrier et al. (1993), p168, ISBN:978-3-540-56670-0\n"
11317 <<
" Söderlind (2002), (2.7) on p5, DOI:10.1023/A:1021160023092\n";
11320 <<
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n"
11337 oomph_info <<
"Tried to increase dt by the ratio "
11340 <<
"). Attempting to increase by the maximum ratio instead."
11353 <<
"Warning: Adaptation of timestep to ensure satisfaction\n"
11354 <<
" of error bounds during adaptive timestepping\n"
11355 <<
" would lower dt below \n"
11356 <<
" Problem::Minimum_dt_but_still_proceed="
11358 <<
" ---> We're continuing with present timestep.\n"
11369 <<
" which is less than the minimum scaling factor "
11371 oomph_info <<
"TIMESTEP REJECTED" << std::endl;
11380 <<
" which is above the maximum (" <<
Maximum_dt
11381 <<
"). I increased it to the maximum value instead.";
11386 std::ostringstream
err;
11388 <<
" which is less than the minimum dt (" <<
Minimum_dt <<
")."
11410 for (
unsigned i = 0;
i <
ni;
i++)
11415#ifdef OOMPH_HAS_MPI
11430 while (reject_timestep);
11457 const double& epsilon,
11477 oomph_info <<
"No spatial refinement allowed; max_adapt=0\n";
11492#ifdef OOMPH_HAS_MPI
11514 oomph_info <<
"Mesh was adapted but re-solve has been suppressed."
11520 <<
"Mesh was adapted --> we'll re-solve for current timestep."
11529 bool shift =
false;
11538 oomph_info <<
"Re-assigning initial condition at time="
11561 oomph_info <<
"Mesh wasn't adapted --> we'll accept spatial refinement."
11661 if (n_sub_mesh == 0)
11688#ifdef OOMPH_HAS_MPI
11692 warning_stream <<
"This has not been comprehensively tested for "
11693 "distributed problems.\n"
11694 <<
"I'm sure that I need to worry about external halo and "
11695 "external elements."
11745 std::string
err =
"Prediction by explicit step only works for "
11746 "problems with a simple time";
11747 err +=
"stepper. I think implementing anything more general will";
11748 err +=
"require a rewrite of explicit time steppers. - David";
11769 std::string
err =
"Requested predictions by explicit step but explicit";
11770 err +=
" predictor pt is null.";
11778 throw OomphLibError(
"Problem has explicit time stepper other than "
11779 "predictor, not sure how to handle this yet ??ds",
11811 using namespace StringConversion;
11812 std::string
err =
"Predictor landed at the wrong time!";
11847#ifdef OOMPH_HAS_MPI
11850 throw OomphLibError(
"Not yet implemented for distributed problems",
11861 std::string
err =
"Not implemented for multiple time steppers";
11871 for (
unsigned i = 0;
i <
ndof();
i++)
11899 elem_pt->enable_mass_matrix_reuse();
11924 elem_pt->disable_mass_matrix_reuse();
11952 oomph_info <<
"Copying an unsteady problem." << std::endl;
11958 for (
unsigned i = 0;
i <
n_dt;
i++)
11978 if (nmesh == 0) nmesh = 1;
11979 for (
unsigned m = 0;
m < nmesh;
m++)
11988 std::ostringstream error_message;
11989 error_message <<
"Number of nodes in copy " <<
n_node
11990 <<
" not equal to the number in the original "
11999 for (
unsigned long i = 0;
i <
n_node;
i++)
12028 std::ostringstream error_message;
12029 error_message <<
"Number of global data in copy " <<
n_global
12030 <<
" not equal to the number in the original "
12047 for (
unsigned m = 0;
m < nmesh;
m++)
12062 std::ostringstream error_message;
12063 error_message <<
"Number of internal data in copy " <<
n_internal
12064 <<
" not equal to the number in the original "
12091 <<
"This function must be overloaded in your specific problem, and must\n"
12092 <<
"create an exact copy of your problem. Usually this will be achieved\n"
12093 <<
"by a call to the constructor with exactly the same arguments as "
12110 dump_file << std::max(
unsigned(1),
n_mesh) <<
" # number of (sub)meshes "
12122 <<
" # uniform refinement when pruned " << std::endl;
12126 dump_file << 0 <<
" # (fake) uniform refinement when pruned "
12129 dump_file << 9999 <<
" # test flag for end of sub-meshes " << std::endl;
12143 <<
" # uniform refinement when pruned " << std::endl;
12147 dump_file << 0 <<
" # (fake) uniform refinement when pruned "
12151 dump_file << 9999 <<
" # test flag for end of sub-meshes " << std::endl;
12154#ifdef OOMPH_HAS_MPI
12162 for (
unsigned e = 0;
e <
n;
e++)
12199 dump_file <<
n <<
" # Number of base elements; partitioning follows.\n";
12200 for (
unsigned e = 0;
e <
n;
e++)
12204 dump_file <<
"8888 #test flag for end of base element distribution\n";
12218#ifdef OOMPH_HAS_TRIANGLE_LIB
12223#ifdef OOMPH_HAS_MPI
12226 if (
mmesh_pt->is_mesh_distributed())
12250#ifdef OOMPH_HAS_TRIANGLE_LIB
12256#ifdef OOMPH_HAS_MPI
12259 if (
mmesh_pt->is_mesh_distributed())
12286 dump_file <<
n_dt <<
" # Number of timesteps " << std::endl;
12287 for (
unsigned i = 0;
i <
n_dt;
i++)
12296 dump_file <<
"0.0 # Dummy time from steady run " << std::endl;
12298 dump_file <<
"0 # Dummy number of timesteps from steady run" << std::endl;
12303 if (nmesh == 0) nmesh = 1;
12304 for (
unsigned m = 0;
m < nmesh;
m++)
12337 warn_message <<
"Restart file isn't open -- I'm assuming that this is\n";
12338 warn_message <<
"because we're restarting on a larger number of\n";
12339 warn_message <<
"processor than were in use when the restart data was \n";
12366 std::ostringstream error_message;
12368 <<
"Number of sub-meshes specified in restart file, "
12369 <<
n_submesh_read <<
" doesn't \n match the my number of sub-meshes,"
12371 <<
"Make sure all sub-meshes have been added to the global mesh\n"
12372 <<
"when calling the Problem::dump() function.\n";
12386#ifdef OOMPH_HAS_MPI
12411 error_stream <<
"Nonzero uniform-refinement-when-pruned specified\n"
12412 <<
"even though mesh is not tree-based. Odd. May want\n"
12413 <<
"to check this carefully before disabling this \n"
12414 <<
"warning/error -- most likely if/when we start to\n"
12415 <<
"prune unstructured meshes [though I can't see why\n"
12416 <<
"we would want to do this, given that they are \n"
12417 <<
"currently totally re-generated...]\n";
12433#ifdef OOMPH_HAS_MPI
12465#ifdef OOMPH_HAS_MPI
12476#ifdef OOMPH_HAS_MPI
12517 std::ostringstream error_message;
12518 error_message <<
"Number of uniform refinements was not consistent \n"
12519 <<
"for following meshes during restart on processor \n"
12520 <<
"on which restart file could be opened:\n";
12527 error_message <<
"Sub-mesh: " <<
i <<
"; local nrefinement: "
12529 <<
"; global/synced nrefinement: "
12558 std::ostringstream error_message;
12560 <<
"Error in reading restart data: Uniform refinement when pruned \n"
12561 <<
"flags should be followed by 9999.\n";
12575#ifdef OOMPH_HAS_MPI
12631 <<
"The number of base elements in the mesh is 0,\n"
12632 <<
" but the restart file indicates that there are "
12634 <<
"This could be because the restart file was \n"
12635 <<
"generated by using code without MPI.\n"
12637 <<
"The best fix is to include two additional lines\n"
12638 <<
"in the restart file: \n\n"
12639 <<
"0 # Number of base elements; partitioning follows.\n"
12640 <<
"8888 # Test flag for end of base element distribution\n"
12642 <<
"These lines go after the flag 9999 that indicates\n"
12643 <<
"the end of the submesh information.\n"
12645 <<
"The file will now continue to be read assuming that\n"
12646 <<
"the base element information is not present.\n"
12647 <<
"If you get strange results then please look carefully\n"
12648 <<
"at the restart file. The safest thing to do is to \n"
12649 <<
"ensure that the restart file was generated by code\n"
12650 <<
"compiled and run with the same parallel options.\n";
12662 std::ostringstream error_message;
12664 <<
" base elements \n"
12665 <<
"though we only have " <<
nbase
12666 <<
" base elements in mesh.\n";
12679 for (
unsigned e = 0;
e <
nbase;
e++)
12703 std::ostringstream error_message;
12705 <<
"Error in reading restart data: Target proc for base elements \n"
12706 <<
"should be followed by 8888.\n";
12727 for (
unsigned e = 0;
e <
nel;
e++)
12756 "el_number_in_base_mesh_plus_one=0 for bulk",
12802 oomph_info <<
"Doing load balancing after pruning\n";
12808 oomph_info <<
"Done load balancing after pruning\n";
12812 oomph_info <<
"No need for load balancing after pruning\n";
12860#ifdef OOMPH_HAS_TRIANGLE_LIB
12865#ifdef OOMPH_HAS_MPI
12868 if (
mmesh_pt->is_mesh_distributed())
12879#ifdef OOMPH_HAS_MPI
12883 if (
mmesh_pt->is_mesh_distributed())
12885 mmesh_pt->reestablish_distribution_info_for_restart(
12916#ifdef OOMPH_HAS_TRIANGLE_LIB
12922#ifdef OOMPH_HAS_MPI
12925 if (
mmesh_pt->is_mesh_distributed())
12937#ifdef OOMPH_HAS_MPI
12941 if (
mmesh_pt->is_mesh_distributed())
12943 mmesh_pt->reestablish_distribution_info_for_restart(
12978 <<
"I've just read in some unstructured meshes and have, in\n"
12979 <<
"the process, totally re-generated their nodes and elements.\n"
12980 <<
"This may create dangling pointers that still point to the\n"
12981 <<
"old nodes and elements, e.g. because FaceElements were\n"
12982 <<
"attached to these meshes or pointers to nodes and elements\n"
12983 <<
"were stored somewhere. FaceElements should therefore be\n"
12984 <<
"removed before reading in these meshes, using an overloaded\n"
12985 <<
"version of the function\n\n"
12986 <<
" Problem::actions_before_read_unstructured_meshes()\n\n"
12987 <<
"and then re-attached using an overloaded version of\n\n"
12988 <<
" Problem::actions_after_read_unstructured_meshes().\n\n"
12989 <<
"The required content of these functions is likely to be "
12991 <<
"to the Problem::actions_before_adapt() and \n"
12992 <<
"Problem::actions_after_adapt() that would be required in\n"
12993 <<
"a spatially adaptive computation. If these functions already\n"
12994 <<
"exist and perform the required actions, the \n"
12995 <<
"actions_before/after_read_unstructured_meshes() functions\n"
12996 <<
"can remain empty because the former are called automatically.\n"
12997 <<
"In this case, this warning my be suppressed by setting the\n"
12998 <<
"public boolean\n\n"
13000 "Problem::Suppress_warning_about_actions_before_read_"
13001 "unstructured_meshes\n\n"
13002 <<
"to true." << std::endl;
13011 oomph_info <<
"\nNumber of equations in Problem::read(): "
13019#ifdef OOMPH_HAS_MPI
13026 oomph_info <<
"Restart file exists" << std::endl;
13027#ifdef OOMPH_HAS_MPI
13074 oomph_info <<
"Restart file does not exist" << std::endl;
13075#ifdef OOMPH_HAS_MPI
13087#ifdef OOMPH_HAS_MPI
13103#ifdef OOMPH_HAS_MPI
13109 std::ostringstream error_message;
13110 error_message <<
"Synchronisation of temporal restart data \n"
13111 <<
"required even though Problem hasn't been distributed "
13172 std::ostringstream error_message;
13174 <<
"Synchronisation of temporal restart data \n"
13175 <<
"required even though we don't have mpi support -- very odd!\n";
13204 if (nmesh == 0) nmesh = 1;
13205 for (
unsigned m = 0;
m < nmesh;
m++)
13253#ifdef OOMPH_HAS_TRIANGLE_LIB
13263 mmesh_pt->update_polyline_representation_from_restart();
13288 std::ostringstream error_message;
13289 error_message <<
"The number of global data " <<
Nglobal
13290 <<
" is not equal to that specified in the input file "
13369 <<
"\n ERROR: Failed Mesh::self_test() for single mesh in problem"
13381 oomph_info <<
"\n ERROR: Failed Mesh::self_test() for mesh imesh"
13382 <<
imesh << std::endl;
13396 <<
"\n ERROR: Failed Data::self_test() for global data iglobal"
13402#ifdef OOMPH_HAS_MPI
13434 const unsigned& bifurcation_type,
13446 double omega = 0.0;
13448 if (bifurcation_type == 3)
13454 double sigma = 0.0;
13455 if (bifurcation_type == 2)
13457 sigma = this->
dof(this->
ndof() - 1);
13471 for (
unsigned c = 0; c <
n_copies; c++)
13492 if (
mmesh_pt->is_adaptation_enabled())
13499 mmesh_pt->refine_base_mesh_as_in_reference_mesh(
13505 <<
"Info/Warning: Mesh in orginal problem is not refineable."
13511 oomph_info <<
"Info/Warning: Mesh adaptation is disabled in copy."
13517 oomph_info <<
"Info/Warning: Mesh cannot be adapted in copy."
13532 if (
mmesh_pt->is_adaptation_enabled())
13539 mmesh_pt->refine_base_mesh_as_in_reference_mesh(
13544 oomph_info <<
"Info/Warning: Mesh in orginal problem is not "
13552 <<
"Info/Warning: Mesh adaptation is disabled in copy."
13558 oomph_info <<
"Info/Warning: Mesh cannot be adapted in copy."
13577 for (
unsigned c = 0; c <
n_copies; c++)
13585 error_stream <<
"Number of unknowns in the problem copy " << c <<
" "
13586 <<
"not equal to number in the original:\n"
13587 << this->
ndof() <<
" (original) "
13602 if (bifurcation_type == 2)
13605 ->symmetrise_eigenfunction_for_adaptive_pitchfork_tracking();
13613 for (
unsigned c = 0; c <
n_copies; c++)
13625 error_stream <<
"Problems do not have the same number of meshes\n"
13643 <<
" does not have the same number of elements in the "
13667 for (
unsigned c = 0; c <
n_copies; c++)
13673 if (bifurcation_type == 2)
13676 ->symmetrise_eigenfunction_for_adaptive_pitchfork_tracking();
13680 for (
unsigned c = 0; c <
n_copies; c++)
13687 switch (bifurcation_type)
13698 this->
dof(this->
ndof() - 1) = sigma;
13709 error_stream <<
"Bifurcation type " << bifurcation_type
13711 <<
"1: Fold, 2: Pitchfork, 3: Hopf\n";
13768 if (bifurcation_type != 0)
13780 for (
unsigned c = 0; c < 2; c++)
13803 if (
mmesh_pt->is_adaptation_enabled())
13813 <<
"Info/Warning: Adaptive Continuation is broken in "
13814 <<
"SolidElement" << std::endl;
13816 mmesh_pt->refine_base_mesh_as_in_reference_mesh(
13821 oomph_info <<
"Info/Warning: Mesh in orginal problem is not "
13829 <<
"Info/Warning: Mesh adaptation is disabled in copy."
13843 <<
"Info/Warning: Adaptive Continuation is broken in "
13844 <<
"SolidElement" << std::endl;
13850 std::ofstream
tri_dump(
"triangle_mesh.dmp");
13853 std::ifstream
tri_read(
"triangle_mesh.dmp");
13867 for (
unsigned i = 0;
i <
n_dim; ++
i)
13876 <<
"Info/warning: Original Mesh is not TriangleBased\n"
13877 <<
"... but the copy is!" << std::endl;
13882 oomph_info <<
"Info/Warning: Mesh cannot be adapted in copy."
13897 if (
mmesh_pt->is_adaptation_enabled())
13907 <<
"Info/Warning: Adaptive Continuation is broken in "
13908 <<
"SolidElement" << std::endl;
13911 mmesh_pt->refine_base_mesh_as_in_reference_mesh(
13916 oomph_info <<
"Info/Warning: Mesh in orginal problem is "
13924 <<
"Info/Warning: Mesh adaptation is disabled in copy."
13938 <<
"Info/Warning: Adaptive Continuation is broken in "
13939 <<
"SolidElement" << std::endl;
13945 std::ofstream
tri_dump(
"triangle_mesh.dmp");
13948 std::ifstream
tri_read(
"triangle_mesh.dmp");
13961 for (
unsigned i = 0;
i <
n_dim; ++
i)
13970 <<
"Info/warning: Original Mesh is not TriangleBased\n"
13971 <<
"... but the copy is!" << std::endl;
13976 oomph_info <<
"Info/Warning: Mesh cannot be adapted in copy."
14003 error_stream <<
"Number of unknowns in the problem copy " << c <<
" "
14004 <<
"not equal to number in the original:\n"
14005 << this->
ndof() <<
" (original) "
14056 oomph_info <<
"Adapting problem:" << std::endl;
14057 oomph_info <<
"=================" << std::endl;
14068 double t_end = 0.0;
14092 if (
mmesh_pt->is_adaptation_enabled())
14098 mmesh_pt->spatial_error_estimator_pt();
14123 mmesh_pt->max_error() = std::fabs(*std::max_element(
14126 mmesh_pt->min_error() = std::fabs(*std::min_element(
14130 <<
mmesh_pt->min_error() << std::endl
14152 oomph_info <<
"Time for complete mesh adaptation "
14153 <<
"(but excluding comp of error estimate): "
14160 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14166 oomph_info <<
"Info/Warning: Mesh cannot be adapted" << std::endl;
14184 mmesh_pt->spatial_error_estimator_pt();
14195 if (
mmesh_pt->is_adaptation_enabled())
14225 <<
mmesh_pt->min_error() << std::endl;
14247 oomph_info <<
"Time for complete mesh adaptation "
14248 <<
"(but excluding comp of error estimate): "
14255 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14261 oomph_info <<
"Info/Warning: Mesh cannot be adapted." << std::endl;
14274 oomph_info <<
"Total time for actual adaptation "
14275 <<
"(all meshes; incl error estimates): " <<
t_end -
t_start
14291 oomph_info <<
"About to start re-assigning eqn numbers "
14292 <<
"with Problem::assign_eqn_numbers() at end of "
14293 <<
"Problem::adapt().\n";
14304 oomph_info <<
"Time for re-assigning eqn numbers with "
14305 <<
"Problem::assign_eqn_numbers() at end of Problem::adapt(): "
14333 if (bifurcation_type != 0)
14341 oomph_info <<
"p-adapting problem:" << std::endl;
14342 oomph_info <<
"===================" << std::endl;
14353 double t_end = 0.0;
14377 if (
mmesh_pt->is_p_adaptation_enabled())
14383 mmesh_pt->spatial_error_estimator_pt();
14408 mmesh_pt->max_error() = std::fabs(*std::max_element(
14411 mmesh_pt->min_error() = std::fabs(*std::min_element(
14415 <<
mmesh_pt->min_error() << std::endl
14437 oomph_info <<
"Time for complete mesh adaptation "
14438 <<
"(but excluding comp of error estimate): "
14445 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14451 oomph_info <<
"Info/Warning: Mesh cannot be adapted" << std::endl;
14469 mmesh_pt->spatial_error_estimator_pt();
14480 if (
mmesh_pt->is_p_adaptation_enabled())
14510 <<
mmesh_pt->min_error() << std::endl;
14532 oomph_info <<
"Time for complete mesh adaptation "
14533 <<
"(but excluding comp of error estimate): "
14540 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14546 oomph_info <<
"Info/Warning: Mesh cannot be adapted." << std::endl;
14559 oomph_info <<
"Total time for actual adaptation "
14560 <<
"(all meshes; incl error estimates): " <<
t_end -
t_start
14576 oomph_info <<
"About to start re-assigning eqn numbers "
14577 <<
"with Problem::assign_eqn_numbers() at end of "
14578 <<
"Problem::adapt().\n";
14589 oomph_info <<
"Time for re-assigning eqn numbers with "
14590 <<
"Problem::assign_eqn_numbers() at end of Problem::adapt(): "
14612 oomph_info <<
"Adapting problem:" << std::endl;
14613 oomph_info <<
"=================" << std::endl;
14633 if (
mmesh_pt->is_adaptation_enabled())
14644 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14650 oomph_info <<
"Info/Warning: Mesh cannot be adapted" << std::endl;
14665 if (
mmesh_pt->is_adaptation_enabled())
14676 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14682 oomph_info <<
"Info/Warning: Mesh cannot be adapted." << std::endl;
14721 if (
mmesh_pt->is_adaptation_enabled())
14725 mmesh_pt->spatial_error_estimator_pt();
14763 <<
mmesh_pt->min_error() << std::endl;
14767 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14773 oomph_info <<
"Info/Warning: Mesh cannot be adapted" << std::endl;
14793 mmesh_pt->spatial_error_estimator_pt();
14804 if (
mmesh_pt->is_adaptation_enabled())
14832 <<
mmesh_pt->min_error() << std::endl;
14836 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14842 oomph_info <<
"Info/Warning: Mesh cannot be adapted." << std::endl;
14857 if (bifurcation_type != 0)
14877 mmesh_pt->spatial_error_estimator_pt();
14901 mmesh_pt->max_error() = std::fabs(*std::max_element(
14904 mmesh_pt->min_error() = std::fabs(*std::min_element(
14908 <<
mmesh_pt->min_error() << std::endl;
14925 mmesh_pt->spatial_error_estimator_pt();
14964 <<
mmesh_pt->min_error() << std::endl;
14995 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15001 std::ostringstream error_message;
15002 error_message <<
"Problem::refine_selected_elements(...) only works for\n"
15003 <<
"multiple-mesh problems if you specify the mesh\n"
15004 <<
"number in the function argument before the Vector,\n"
15005 <<
"or a Vector of Vectors for each submesh.\n"
15041 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15047 std::ostringstream error_message;
15048 error_message <<
"Problem::refine_selected_elements(...) only works for\n"
15049 <<
"multiple-mesh problems if you specify the mesh\n"
15050 <<
"number in the function argument before the Vector,\n"
15051 <<
"or a Vector of Vectors for each submesh.\n"
15078 std::ostringstream error_message;
15079 error_message <<
"Problem only has " <<
n_mesh
15080 <<
" submeshes. Cannot refine submesh " <<
i_mesh
15094 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15126 std::ostringstream error_message;
15127 error_message <<
"Problem only has " <<
n_mesh
15128 <<
" submeshes. Cannot refine submesh " <<
i_mesh
15142 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15181 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15218 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15256 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15262 std::ostringstream error_message;
15264 <<
"Problem::p_refine_selected_elements(...) only works for\n"
15265 <<
"multiple-mesh problems if you specify the mesh\n"
15266 <<
"number in the function argument before the Vector,\n"
15267 <<
"or a Vector of Vectors for each submesh.\n"
15303 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15309 std::ostringstream error_message;
15311 <<
"Problem::p_refine_selected_elements(...) only works for\n"
15312 <<
"multiple-mesh problems if you specify the mesh\n"
15313 <<
"number in the function argument before the Vector,\n"
15314 <<
"or a Vector of Vectors for each submesh.\n"
15335 "p-refinement for multiple submeshes has not yet been tested.",
15336 "Problem::p_refine_selected_elements()",
15346 std::ostringstream error_message;
15347 error_message <<
"Problem only has " <<
n_mesh
15348 <<
" submeshes. Cannot p-refine submesh " <<
i_mesh
15362 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15388 "p-refinement for multiple submeshes has not yet been tested.",
15389 "Problem::p_refine_selected_elements()",
15399 std::ostringstream error_message;
15400 error_message <<
"Problem only has " <<
n_mesh
15401 <<
" submeshes. Cannot p-refine submesh " <<
i_mesh
15415 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15440 "p-refinement for multiple submeshes has not yet been tested.",
15441 "Problem::p_refine_selected_elements()",
15459 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15482 "p-refinement for multiple submeshes has not yet been tested.",
15483 "Problem::p_refine_selected_elements()",
15501 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15535 double t_end = 0.0;
15540 <<
"Time for actions before adapt in Problem::refine_uniformly_aux(): "
15556 for (
unsigned i = 0;
i <
nref;
i++)
15558 mmesh_pt->refine_uniformly(doc_info);
15563 oomph_info <<
"Info/Warning: Mesh cannot be refined uniformly "
15578 for (
unsigned i = 0;
i <
nref;
i++)
15580 mmesh_pt->refine_uniformly(doc_info);
15596 oomph_info <<
"Time for mesh-level mesh refinement in "
15597 <<
"Problem::refine_uniformly_aux(): " <<
t_end -
t_start
15610 <<
"Time for actions after adapt Problem::refine_uniformly_aux(): "
15616#ifdef OOMPH_HAS_MPI
15629 oomph_info <<
"Time for Problem::prune_halo_elements_and_nodes() in "
15630 <<
"Problem::refine_uniformly_aux(): " <<
t_end -
t_start
15638 std::ostringstream error_message;
15640 <<
"Requested pruning in serial build. Ignoring the request.\n";
15642 "Problem::refine_uniformly_aux()",
15649 <<
"Number of equations after Problem::refine_uniformly_aux(): "
15655 oomph_info <<
"Time for Problem::assign_eqn_numbers() in "
15656 <<
"Problem::refine_uniformly_aux(): " <<
t_end -
t_start
15682 double t_end = 0.0;
15686 oomph_info <<
"Time for actions before adapt in "
15687 "Problem::p_refine_uniformly_aux(): "
15703 for (
unsigned i = 0;
i <
nref;
i++)
15705 mmesh_pt->p_refine_uniformly(doc_info);
15710 oomph_info <<
"Info/Warning: Mesh cannot be p-refined uniformly "
15718 "p-refinement for multiple submeshes has not yet been tested.",
15719 "Problem::p_refine_uniformly_aux()",
15730 for (
unsigned i = 0;
i <
nref;
i++)
15732 mmesh_pt->p_refine_uniformly(doc_info);
15748 oomph_info <<
"Time for mesh-level mesh refinement in "
15749 <<
"Problem::p_refine_uniformly_aux(): " <<
t_end -
t_start
15762 <<
"Time for actions after adapt Problem::p_refine_uniformly_aux(): "
15768#ifdef OOMPH_HAS_MPI
15781 oomph_info <<
"Time for Problem::prune_halo_elements_and_nodes() in "
15782 <<
"Problem::p_refine_uniformly_aux(): " <<
t_end -
t_start
15790 std::ostringstream error_message;
15792 <<
"Requested pruning in serial build. Ignoring the request.\n";
15794 "Problem::p_refine_uniformly_aux()",
15801 <<
"Number of equations after Problem::p_refine_uniformly_aux(): "
15807 oomph_info <<
"Time for Problem::assign_eqn_numbers() in "
15808 <<
"Problem::p_refine_uniformly_aux(): " <<
t_end -
t_start
15826 std::ostringstream error_message;
15827 error_message <<
"imesh " <<
i_mesh
15828 <<
" is greater than the number of sub meshes "
15840 mmesh_pt->refine_uniformly(doc_info);
15844 oomph_info <<
"Info/Warning: Mesh cannot be refined uniformly "
15870 std::ostringstream error_message;
15871 error_message <<
"imesh " <<
i_mesh
15872 <<
" is greater than the number of sub meshes "
15884 mmesh_pt->p_refine_uniformly(doc_info);
15888 oomph_info <<
"Info/Warning: Mesh cannot be refined uniformly "
15931 oomph_info <<
"Info/Warning: Mesh cannot be unrefined uniformly "
15991 std::ostringstream error_message;
15992 error_message <<
"imesh " <<
i_mesh
15993 <<
" is greater than the number of sub meshes "
16009 oomph_info <<
"Info/Warning: Mesh cannot be unrefined uniformly "
16052 mmesh_pt->p_unrefine_uniformly(doc_info);
16056 oomph_info <<
"Info/Warning: Mesh cannot be p-unrefined uniformly "
16064 throw OomphLibError(
"This functionality has not yet been tested.",
16074 mmesh_pt->p_unrefine_uniformly(doc_info);
16105 std::ostringstream error_message;
16106 error_message <<
"imesh " <<
i_mesh
16107 <<
" is greater than the number of sub meshes "
16119 mmesh_pt->p_unrefine_uniformly(doc_info);
16123 oomph_info <<
"Info/Warning: Mesh cannot be p-unrefined uniformly "
16163 <<
"\n\n===========================================================\n";
16164 oomph_info <<
" ******** WARNING *********** \n";
16166 <<
"===========================================================\n";
16167 oomph_info <<
"Problem::unsteady_newton_solve() called with "
16171 oomph_info <<
"This doesn't make sense (shifting does have to be done"
16173 oomph_info <<
"since we're constantly re-assigning the initial conditions"
16176 <<
"\n===========================================================\n\n";
16199#ifdef OOMPH_HAS_MPI
16224 <<
n_unrefined <<
" were unrefined, in total." << std::endl;
16229 oomph_info <<
"\n \n Solution is fully converged in "
16230 <<
"Problem::unsteady_newton_solver() \n \n ";
16245 oomph_info <<
"Re-setting initial condition " << std::endl;
16280 <<
"----------------------------------------------------------"
16282 <<
"Reached max. number of adaptations in \n"
16283 <<
"Problem::unsteady_newton_solver().\n"
16284 <<
"----------------------------------------------------------"
16319#ifdef OOMPH_HAS_MPI
16345#ifdef OOMPH_HAS_MPI
16346 <<
", in total (over all processors).\n";
16355 oomph_info <<
"\n \n Solution is fully converged in "
16356 <<
"Problem::newton_solver(). \n \n ";
16378 <<
"USER-DEFINED ERROR IN NEWTON SOLVER " << std::endl;
16383 <<
") REACHED WITHOUT CONVERGENCE " << std::endl;
16395 error_stream <<
"Error occured in adaptive Newton solver. "
16413 <<
"----------------------------------------------------------"
16415 <<
"Reached max. number of adaptations in \n"
16416 <<
"Problem::newton_solver().\n"
16417 <<
"----------------------------------------------------------"
16449#ifdef OOMPH_HAS_MPI
16492 oomph_info <<
"Checking halo schemes on single mesh" << std::endl;
16493 doc_info.
label() =
"_one_and_only_mesh_";
16503 std::stringstream
tmp;
16601 for (
unsigned n = 0;
n <
n_nod;
n++)
16727 for (
unsigned n = 0;
n <
n_nod;
n++)
16838 double t_end = 0.0;
16864 for (
unsigned e = 0;
e <
nelem;
e++)
16889 for (
unsigned j = 0;
j <
nnod;
j++)
16950 oomph_info <<
"Time for copy_haloed_eqn_numbers_helper for halos: "
16964 <<
"Time for copy_haloed_eqn_numbers_helper for external halos: "
16973 if (assign_local_eqn_numbers)
16995 oomph_info <<
"Time for assign_local_eqn_numbers in sync: "
17088 for (
unsigned n = 0;
n <
n_nod;
n++)
17112 ->add_eqn_numbers_to_vector(
send_data);
17215 for (
unsigned n = 0;
n <
n_nod;
n++)
17284 warn_message <<
"WARNING: You've tried to load balance a problem over\n"
17285 <<
"only one processor: ignoring your request.\n";
17287 "Problem::load_balance()",
17299 error_stream <<
"You have called Problem::load_balance()\n"
17300 <<
"on a non-distributed problem. This doesn't\n"
17301 <<
"make sense -- go distribute your problem first."
17409 std::map<GeneralisedElement*, unsigned>
17496 <<
"The number of nonhalo elements (" <<
count_td
17497 <<
") found in (all)\n"
17498 <<
"the (sub)-mesh(es) is different from the number of target "
17501 <<
"Please ensure that you called the rebuild_global_mesh() method "
17502 <<
"after the\npossible deletion of FaceElements in "
17503 <<
"actions_before_distribute()!!!\n\n";
17505 "Problem::load_balance()",
17671 oomph_info <<
"CPU for partition calculation for roots: "
17694 ref_mesh_pt->uniform_refinement_level_when_pruned();
17720 ref_mesh_pt->uniform_refinement_level_when_pruned();
17750 for (
unsigned i = 0;
i <
n;
i++)
17808 for (
unsigned i = 0;
i <
nref;
i++)
17851 <<
nref <<
" times\n";
17852 for (
unsigned i = 0;
i <
nref;
i++)
17892 unsigned count = 0;
17916 <<
"The number of READ target domains for nonhalo elements\n"
17917 <<
" is (" <<
count <<
"), but the number of target domains for\n"
17921 "Problem::load_balance()",
17941 for (
unsigned e = 0;
e <
n_ele;
e++)
17975 for (
unsigned e = 0;
e <
n_ele;
e++)
17991 <<
") is not the same as the number of\nadded elements ("
17992 <<
counter_base <<
") to the Base_mesh_element_pt data "
17993 <<
"structure!!!\n\n";
17995 "Problem::load_balance()",
18022 std::stringstream
info;
18023 info <<
"================================================\n";
18024 info <<
"INFO: I've come across a FaceElement while \n";
18025 info <<
" load-balancing a problem. \n";
18026 info <<
"================================================\n";
18034 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
18058 error_stream <<
"Distributing one-and-only mesh containing "
18069 oomph_info <<
"Distributing one and only mesh\n"
18070 <<
"------------------------------" << std::endl;
18078 new_domain_for_base_element,
18100 <<
n_mesh <<
" in total\n"
18101 <<
"---------------------------------------------"
18112 submesh_element_partition[
i_mesh],
18135 for (
unsigned e = 0;
e <
n_del;
e++)
18176 for (
unsigned e = 0;
e <
n_del;
e++)
18191 oomph_info <<
"CPU for build and distribution of new mesh(es): "
18225 oomph_info <<
"CPU for refinement of base mesh: "
18271 this->
mesh_pt() = old_mesh_pt[0];
18277 error_stream <<
"The only one mesh in the problem is not an "
18278 "unstructured mesh,\n"
18279 <<
"but the flag 'are_there_unstructures_meshes' ("
18281 <<
") was turned on,\n"
18282 <<
"this is weird. Please check for any condition "
18284 <<
"turned on this flag!!!!\n\n";
18286 "Problem::load_balance()",
18350 oomph_info <<
"CPU for transferring solution to new mesh(es): "
18369 <<
"Total number of elements on this processor after load balance: "
18372 oomph_info <<
"Number of non-halo elements on this processor after "
18382 <<
"Number of dofs in load_balance() has changed from " <<
old_ndof
18383 <<
" to " <<
n_dof <<
"\n"
18384 <<
"Check that you've implemented any necessary "
18385 "actions_before/after\n"
18386 <<
"adapt/distribute functions, e.g. to pin redundant pressure dofs"
18460 for (
unsigned h = 0; h <
nhaloed; h++)
18467 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
18492 unsigned count = 0;
18519 for (
unsigned j = 0;
j <
nhalo;
j++)
18576 std::stringstream error_message;
18578 <<
" doesn't match that actually send: "
18731 for (
unsigned j = 0;
j <
nhalo;
j++)
18814 std::stringstream error_message;
18816 <<
" doesn't match that actually send: "
18854 for (
unsigned jd = 0;
jd <
nd;
jd++)
18874 for (
unsigned j = 0;
j <
n;
j++)
19006 for (
unsigned j = 0;
j <
n;
j++)
19034 const int n_proc = comm_pt->nproc();
19107 for (
unsigned b = 0; b <
nbatch; b++)
19140 std::ostringstream error_message;
19142 <<
"Number of leaves: " <<
n_leaf <<
" "
19143 <<
" doesn't match number of elements sent in batch: " <<
nel
19165 std::ostringstream error_message;
19167 <<
"Non-refineable root element should only be associated with"
19168 <<
" one element but nel=" <<
nel <<
"\n";
19178 for (
unsigned e = 0;
e <
nel;
e++)
19189 for (
unsigned j = 0;
j <
nnod;
j++)
19207 std::ostringstream error_message;
19209 <<
"Node has more values, namely " <<
nod_pt->nvalue()
19210 <<
", than we're about to receive, namely " <<
nval
19211 <<
". Something's wrong!\n";
19234 std::ostringstream error_message;
19235 error_message <<
"Local node is boundary node but "
19236 "information sent is\n"
19237 <<
"for non-boundary node\n";
19254 ->index_of_first_value_assigned_by_face_element_pt() ==
19258 ->index_of_first_value_assigned_by_face_element_pt() =
19259 new std::map<unsigned, unsigned>;
19264 std::map<unsigned, unsigned>*
map_pt =
19266 ->index_of_first_value_assigned_by_face_element_pt();
19289 std::ostringstream error_message;
19290 error_message <<
"Local node is not a boundary node but "
19292 <<
"sent is for boundary node.\n";
19368 const int n_proc = comm_pt->nproc();
19389 unsigned max_refinement_level = 0;
19410 std::map<unsigned, Vector<unsigned>>
19459 for (
unsigned e = 0;
e <
nele;
e++)
19491 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
19528 "Base_mesh_element_number_plus_one[...]=0",
19550 root_el_pt->tree_pt()->stick_leaves_into_vector(
19574 if (level > max_refinement_level)
19576 max_refinement_level = level;
19590 std::ostringstream error_message;
19592 <<
"All elements associated with same root must have "
19593 <<
"same target. during load balancing\n";
19610 std::ostringstream error_message;
19614 <<
" target domains for local non-halo elelemts. \n "
19615 <<
"Very Odd -- we do (now) strip out the information for elements\n"
19616 <<
"that are removed in actions_before_distribute()...\n";
19639 comm_pt->mpi_comm());
19655 comm_pt->mpi_comm());
19668 comm_pt->mpi_comm());
19679 std::ostringstream error_message;
19680 error_message <<
"Old domain for base element " <<
j <<
": "
19682 <<
"or its incarnation as refineable el: "
19685 <<
" which is of type "
19688 <<
"appear to have been assigned by any processor\n";
19698 std::ostringstream error_message;
19699 error_message <<
"New domain for base element " <<
j
19700 <<
"which is of type "
19703 <<
"appear to have been assigned by any processor\n";
19744 for (
unsigned b = 0; b <
nbatch; b++)
19760 for (
unsigned e = 0;
e <
nel;
e++)
19771 for (
unsigned j = 0;
j <
nnod;
j++)
19778 unsigned nt =
nod_pt->ntstorage();
19784 for (
unsigned t = 0;
t <
nt;
t++)
19792 for (
unsigned i = 0;
i <
n_dim;
i++)
19831 std::map<unsigned, unsigned>*
map_pt =
19832 bnod_pt->index_of_first_value_assigned_by_face_element_pt();
19846 for (std::map<unsigned, unsigned>::iterator
p =
19851 send_data.push_back(
double((*p).first));
19852 send_data.push_back(
double((*p).second));
19877 std::ostringstream error_message;
19880 <<
" doesn't match total number of elements sent in batch: "
19958 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
19983 "Base_mesh_element_number_plus_one[...]=0",
19992 root_el_pt->tree_pt()->stick_all_tree_nodes_into_vector(
20117 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
20232 for (
unsigned e = 0;
e <
nel;
e++)
20309 for (
unsigned e = 0;
e <
nel;
e++)
A class that is used to define the functions used to assemble the elemental contributions to the resi...
virtual unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
virtual void synchronise()
Function that is used to perform any synchronisation required during the solution.
virtual int bifurcation_type() const
Return an unsigned integer to indicate whether the handler is a bifurcation tracking handler....
virtual void get_hessian_vector_products(GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Calculate the product of the Hessian (derivative of Jacobian with respect to all variables) an eigenv...
virtual double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
virtual void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
virtual void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt.
virtual unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Return the global equation number of the local unknown ieqn_local in elem_pt.
virtual void get_all_vectors_and_matrices(GeneralisedElement *const &elem_pt, Vector< Vector< double > > &vec, Vector< DenseMatrix< double > > &matrix)
Calculate all desired vectors and matrices provided by the element elem_pt.
virtual void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable" for elem_pt.
A custom linear solver class that is used to solve a block-factorised version of the Hopf bifurcation...
A class that contains the information required by Nodes that are located on Mesh boundaries....
A class for compressed column matrices that store doubles.
void build_without_copy(T *value, int *row_index, int *column_start, const unsigned long &nnz, const unsigned long &n, const unsigned long &m)
Function to build matrix from pointers to arrays which hold the column starts, row indices and non-ze...
A class for compressed row matrices. This is a distributable object.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the matrix are redistributed to match the new distribution. In a non-MPI build this m...
void build_without_copy(const unsigned &ncol, const unsigned &nnz, double *value, int *column_index, int *row_start)
keeps the existing distribution and just matrix that is stored without copying the matrix data
void build(const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)
build method: vector of values, vector of column indices, vector of row starts and number of rows and...
A Base class for DGElements.
A class that represents a collection of data; each Data object may contain many different individual ...
void copy(Data *orig_data_pt)
Copy Data values from specified Data object.
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
void initialise(const T &val)
Initialize all values in the matrix to val.
void resize(const unsigned long &n)
Resize to a square nxn matrix; any values already present will be transfered.
bool distribution_built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
Information for documentation of results: Directory and file number to enable output in the form RESL...
std::string & label()
String used (e.g.) for labeling output files.
bool is_doc_enabled() const
Are we documenting?
void disable_doc()
Disable documentation.
std::string directory() const
Output directory.
unsigned & number()
Number used (e.g.) for labeling output files.
A class that stores the halo/haloed entries required when using a DoubleVectorWithHaloEntries....
void setup_halo_dofs(const std::map< unsigned, double * > &halo_data_pt, Vector< double * > &halo_dof_pt)
Function that sets up a vector of pointers to halo data, index using the scheme in Local_index.
===================================================================== An extension of DoubleVector th...
void build_halo_scheme(DoubleVectorHaloScheme *const &halo_scheme_pt)
Construct the halo scheme and storage for the halo data.
void sum_all_halo_and_haloed_values()
Sum all the data, store in the master (haloed) data and then synchronise.
double & global_value(const unsigned &i)
Direct access to global entry.
A vector in the mathematical sense, initially developed for linear algebra type applications....
double max() const
returns the maximum coefficient
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
double * values_pt()
access function to the underlying values
void clear()
wipes the DoubleVector
A class that is used to define the functions used to assemble the elemental contributions to the mass...
virtual void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector_real, Vector< DoubleVector > &eigenvector_imag, const bool &do_adjoint_problem=false)
Solve the real eigenproblem that is assembled by elements in a mesh in a Problem object....
Base class for spatial error estimators.
A class that is used to define the functions used to assemble and invert the mass matrix when taking ...
A Base class for explicit timesteppers.
virtual void timestep(ExplicitTimeSteppableObject *const &object_pt, const double &dt)=0
Pure virtual function that is used to advance time in the object.
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
A general Finite Element class.
void position(const Vector< double > &zeta, Vector< double > &r) const
Return the parametrised position of the FiniteElement in its incarnation as a GeomObject,...
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
unsigned nnode() const
Return the number of nodes.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
virtual void describe_local_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the local dofs of the element[s]. The ostream specifies the output stream to whi...
A Generalised Element class.
bool is_halo() const
Is this element a halo?
void read_internal_eqn_numbers_from_vector(const Vector< long > &vector_of_eqn_numbers, unsigned &index)
Read all equation numbers associated with internal data from the vector starting from index....
unsigned ndof() const
Return the number of equations/dofs in the element.
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
void add_internal_data_values_to_vector(Vector< double > &vector_of_values)
Add all internal data and time history values to the vector in the internal storage order.
unsigned ninternal_data() const
Return the number of internal data objects.
void add_internal_eqn_numbers_to_vector(Vector< long > &vector_of_eqn_numbers)
Add all equation numbers associated with internal data to the vector in the internal storage order.
virtual void assign_local_eqn_numbers(const bool &store_local_dof_pt)
Setup the arrays of local equation numbers for the element. If the optional boolean argument is true,...
void describe_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the dofs of the element. The ostream specifies the output stream to which the de...
virtual void complete_setup_of_dependencies()
Complete the setup of any additional dependencies that the element may have. Empty virtual function t...
void read_internal_data_values_from_vector(const Vector< double > &vector_of_values, unsigned &index)
Read all internal data and time history values from the vector starting from index....
unsigned ndim() const
Access function to # of Eulerian coordinates.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Class that contains data for hanging nodes.
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
unsigned nmaster() const
Return the number of master nodes.
void set_master_node_pt(const unsigned &i, Node *const &master_node_pt, const double &weight)
Set the pointer to the i-th master node and its weight.
A class that is used to assemble the augmented system that defines a Hopf bifurcation....
A class to specify when the error is caused by an inverted element.
Class for the LAPACK QZ eigensolver.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
void build(const OomphCommunicator *const comm_pt, const unsigned &first_row, const unsigned &nrow_local, const unsigned &nrow=0)
Sets the distribution. Takes first_row, nrow_local and nrow as arguments. If nrow is not provided or ...
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
virtual void enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
virtual void enable_computation_of_gradient()
function to enable the computation of the gradient required for the globally convergent Newton method
virtual void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled jacobian and the rhs vector. Solution is returned in...
void get_gradient(DoubleVector &gradient)
function to access the gradient, provided it has been computed
void reset_gradient()
function to reset the size of the gradient before each Newton solve
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
bool is_resolve_enabled() const
Boolean flag indicating if resolves are enabled.
static bool mpi_has_been_initialised()
return true if MPI has been initialised
static OomphCommunicator * communicator_pt()
access to global communicator. This is the oomph-lib equivalent of MPI_COMM_WORLD
bool does_pointer_correspond_to_mesh_data(double *const ¶meter_pt)
Does the double pointer correspond to any mesh data.
void remove_boundary_node(const unsigned &b, Node *const &node_pt)
Remove a node from the boundary b.
GeneralisedElement *& external_halo_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th external halo element in this Mesh whose non-halo counterpart is held on proce...
void flush_element_and_node_storage()
Flush storage for elements and nodes by emptying the vectors that store the pointers to them....
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
void check_halo_schemes(DocInfo &doc_info, double &max_permitted_error_for_halo_check)
Check halo and shared schemes on the mesh.
virtual void set_mesh_level_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Function that can be used to set any additional timestepper data stored at the Mesh (as opposed to no...
void describe_local_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the local dofs of the elements. The ostream specifies the output stream to which...
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
void set_nodal_and_elemental_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with all nodal and elemental data stored in the mesh.
void shift_time_values()
Shift time-dependent data along for next timestep: Deal with nodal Data/positions and the element's i...
void prune_halo_elements_and_nodes(Vector< GeneralisedElement * > &deleted_element_pt, const bool &report_stats=false)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement,...
void calculate_predictions()
Calculate predictions for all Data and positions associated with the mesh, usually used in adaptive t...
void describe_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the dofs of the Mesh. The ostream specifies the output stream to which the descr...
unsigned long nnode() const
Return number of nodes in the mesh.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
void delete_all_external_storage()
Wipe the storage for all externally-based elements.
unsigned nnon_halo_element()
Total number of non-halo elements in this mesh (Costly call computes result on the fly)
void get_all_halo_data(std::map< unsigned, double * > &map_of_halo_data)
Get all the halo data stored in the mesh and add pointers to the data to the map, indexed by global e...
void assign_initial_values_impulsive()
Assign initial values for an impulsive start.
void assign_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers in all elements If the boolean argument is true then also store poi...
virtual void read(std::ifstream &restart_file)
Read solution from restart file.
void set_consistent_pinned_values_for_continuation(ContinuationStorageScheme *const &continuation_stepper_pt)
Set consistent values for pinned data in continuation.
unsigned long assign_global_eqn_numbers(Vector< double * > &Dof_pt)
Assign the global equation numbers in the Data stored at the nodes and also internal element Data....
virtual void distribute(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, Vector< GeneralisedElement * > &deleted_element_pt, DocInfo &doc_info, const bool &report_stats, const bool &overrule_keep_as_halo_element_status)
Distribute the problem and doc; make this virtual to allow overloading for particular meshes where fu...
void null_external_halo_node(const unsigned &p, Node *nod_pt)
Null out specified external halo node (used when deleting duplicates)
unsigned long nelement() const
Return number of elements in the mesh.
void merge_meshes(const Vector< Mesh * > &sub_mesh_pt)
Merge meshes. Note: This simply merges the meshes' elements and nodes (ignoring duplicates; no bounda...
unsigned nexternal_halo_element()
Total number of external halo elements in this Mesh.
virtual void dump(std::ofstream &dump_file, const bool &use_old_ordering=true) const
Dump the data in the mesh into a file for restart.
A class to handle errors in the Newton solver.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
void copy(Node *orig_node_pt)
Copy all nodal data from specified Node object.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
bool is_hanging() const
Test whether the node is geometrically hanging.
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
std::ostream *& stream_pt()
Access function for the stream pointer.
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
A class that is used to assemble the residuals in parallel by overloading the get_all_vectors_and_mat...
A class that is used to define the functions used when assembling the derivatives of the residuals wi...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
virtual void actions_after_implicit_timestep()
Actions that should be performed after each implicit time step. This is needed when one wants to solv...
bool Always_take_one_newton_step
Boolean to indicate whether a Newton step should be taken even if the initial residuals are below the...
bool Jacobian_reuse_is_enabled
Is re-use of Jacobian in Newton iteration enabled? Default: false.
virtual void actions_after_newton_solve()
Any actions that are to be performed after a complete Newton solve, e.g. post processing....
void parallel_sparse_assemble(const LinearAlgebraDistribution *const &dist_pt, Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residuals)
Helper method to assemble CRDoubleMatrices from distributed on multiple processors.
void describe_dofs(std::ostream &out= *(oomph_info.stream_pt())) const
Function to describe the dofs in terms of the global equation number, i.e. what type of value (nodal ...
virtual void sparse_assemble_row_or_column_compressed(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Protected helper function that is used to assemble the Jacobian matrix in the case when the storage i...
double * global_dof_pt(const unsigned &i)
Return a pointer to the dof, indexed by global equation number which may be haloed or stored locally....
bool Store_local_dof_pt_in_elements
Boolean to indicate whether local dof pointers should be stored in the elements.
virtual void actions_before_newton_step()
Any actions that are to be performed before each individual Newton step. Most likely to be used for d...
void remove_duplicate_data(Mesh *const &mesh_pt, bool &actually_removed_some_data)
Private helper function to remove repeated data in external haloed elements in specified mesh....
bool Bifurcation_detection
Boolean to control bifurcation detection via determinant of Jacobian.
void get_flat_packed_refinement_pattern_for_load_balancing(const Vector< unsigned > &old_domain_for_base_element, const Vector< unsigned > &new_domain_for_base_element, const unsigned &max_refinement_level_overall, std::map< unsigned, Vector< unsigned > > &flat_packed_refinement_info_for_root)
Get flat-packed refinement pattern for each root element in current mesh (labeled by unique number of...
void adapt()
Adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on their own error e...
virtual void actions_before_newton_solve()
Any actions that are to be performed before a complete Newton solve (e.g. adjust boundary conditions)...
Vector< unsigned > First_el_for_assembly
First element to be assembled by given processor for non-distributed problem (only kept up to date wh...
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
void refine_uniformly_aux(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info, const bool &prune)
Helper function to do compund refinement of (all) refineable (sub)mesh(es) uniformly as many times as...
void build_global_mesh()
Build the global mesh by combining the all the submeshes. Note: The nodes boundary information refers...
unsigned unrefine_uniformly()
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem. Return 0 for success,...
void assign_initial_values_impulsive()
Initialise data and nodal positions to simulate impulsive start from initial configuration/solution.
double Theta_squared
Value of the scaling parameter required so that the parameter occupies the desired proportion of the ...
virtual void get_eigenproblem_matrices(CRDoubleMatrix &mass_matrix, CRDoubleMatrix &main_matrix, const double &shift=0.0)
Get the matrices required by a eigensolver. If the shift parameter is non-zero the second matrix will...
Vector< double > Dof_derivative
Storage for the derivative of the problem variables wrt arc-length.
double & dof_current(const unsigned &i)
Access function to the current value of the i-th (local) dof at the start of a continuation step.
virtual void sparse_assemble_row_or_column_compressed_with_lists(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
virtual void sparse_assemble_row_or_column_compressed_with_two_arrays(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
friend class BlockHopfLinearSolver
void synchronise_dofs(const bool &do_halos, const bool &do_external_halos)
Synchronise the degrees of freedom by overwriting the haloed values with their non-halo counterparts ...
virtual void actions_before_distribute()
Actions to be performed before a (mesh) distribution.
Distributed_problem_matrix_distribution Dist_problem_matrix_distribution
The distributed matrix distribution method 1 - Automatic - the Problem distribution is employed,...
DoubleVectorHaloScheme * Halo_scheme_pt
Pointer to the halo scheme for any global vectors that have the Dof_distribution.
virtual void actions_after_change_in_global_parameter(double *const ¶meter_pt)
Actions that are to be performed when the global parameter addressed by parameter_pt has been changed...
void p_refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
p-refine (one and only!) mesh by refining the elements identified by their numbers relative to the pr...
void setup_dof_halo_scheme()
Function that is used to setup the halo scheme.
unsigned long ndof() const
Return the number of dofs.
void p_unrefine_uniformly(DocInfo &doc_info)
p-unrefine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem.
virtual void build_mesh()
Function to build the Problem's base mesh; this must be supplied by the user if they wish to use the ...
static bool Suppress_warning_about_actions_before_read_unstructured_meshes
Flag to allow suppression of warning messages re reading in unstructured meshes during restart.
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
bool Use_default_partition_in_load_balance
Flag to use "default partition" during load balance. Should only be set to true when run in validatio...
void bifurcation_adapt_helper(unsigned &n_refined, unsigned &n_unrefined, const unsigned &bifurcation_type, const bool &actually_adapt=true)
A function that is used to adapt a bifurcation-tracking problem, which requires separate interpolatio...
void adapt_based_on_error_estimates(unsigned &n_refined, unsigned &n_unrefined, Vector< Vector< double > > &elemental_error)
Adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on the error estimat...
friend class AugmentedBlockFoldLinearSolver
Vector< double > Elemental_assembly_time
Storage for assembly times (used for load balancing)
double & dof(const unsigned &i)
i-th dof in the problem
bool Jacobian_has_been_computed
Has a Jacobian been computed (and can therefore be re-used if required)? Default: false.
bool Bypass_increase_in_dof_check_during_pruning
Boolean to bypass check of increase in dofs during pruning.
virtual void get_dvaluesdt(DoubleVector &f)
Get the time derivative of all values (using get_inverse_mass_matrix_times_residuals(....
void initialise_dt(const double &dt)
Set all timesteps to the same value, dt, and assign weights for all timesteppers in the problem.
void add_to_dofs(const double &lambda, const DoubleVector &increment_dofs)
Add lambda x incremenet_dofs[l] to the l-th dof.
double Timestep_reduction_factor_after_nonconvergence
What it says: If temporally adaptive Newton solver fails to to converge, reduce timestep by this fact...
Vector< double > Dof_current
Storage for the present values of the variables.
double Desired_proportion_of_arc_length
Proportion of the arc-length to taken by the parameter.
virtual void actions_after_implicit_timestep_and_error_estimation()
Actions that should be performed after each implicit time step. This is needed if your actions_after_...
void disable_mass_matrix_reuse()
Turn off recyling of the mass matrix in explicit timestepping schemes.
EigenSolver * Default_eigen_solver_pt
Pointer to the default eigensolver.
unsigned Nnewton_iter_taken
Actual number of Newton iterations taken during the most recent iteration.
double * bifurcation_parameter_pt() const
Return pointer to the parameter that is used in the bifurcation detection. If we are not tracking a b...
void copy_haloed_eqn_numbers_helper(const bool &do_halos, const bool &do_external_halos)
A private helper function to copy the haloed equation numbers into the halo equation numbers,...
virtual void sparse_assemble_row_or_column_compressed_with_vectors_of_pairs(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
void send_data_to_be_sent_during_load_balancing(Vector< int > &send_n, Vector< double > &send_data, Vector< int > &send_displacement)
Load balance helper routine: Send data to other processors during load balancing.
bool Time_adaptive_newton_crash_on_solve_fail
Bool to specify what to do if a Newton solve fails within a time adaptive solve. Default (false) is t...
Problem()
Constructor: Allocate space for one time stepper and set all pointers to NULL and set defaults for al...
bool is_dparameter_calculated_analytically(double *const ¶meter_pt)
Function to determine whether the parameter derivatives are calculated analytically.
void flush_sub_meshes()
Flush the problem's collection of sub-meshes. Must be followed by call to rebuild_global_mesh().
unsigned Sparse_assembly_method
Flag to determine which sparse assembly method to use By default we use assembly by vectors of pairs.
bool & use_predictor_values_as_initial_guess()
void calculate_predictions()
Calculate predictions.
Vector< unsigned > Last_el_plus_one_for_assembly
Last element (plus one) to be assembled by given processor for non-distributed problem (only kept up ...
virtual void actions_after_read_unstructured_meshes()
Actions that are to be performed before reading in restart data for problems involving unstructured b...
void copy(Problem *orig_problem_pt)
Copy Data values, nodal positions etc from specified problem. Note: This is not a copy constructor....
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
void set_explicit_time_stepper_pt(ExplicitTimeStepper *const &explicit_time_stepper_pt)
Set the explicit timestepper for the problem. The function will automatically create or resize the Ti...
virtual void actions_after_distribute()
Actions to be performed after a (mesh) distribution.
virtual void actions_before_implicit_timestep()
Actions that should be performed before each implicit time step. This is needed when one wants to sol...
LinearSolver * Linear_solver_pt
Pointer to the linear solver for the problem.
bool Doc_time_in_distribute
Protected boolean flag to provide comprehensive timimings during problem distribution....
unsigned Max_newton_iterations
Maximum number of Newton iterations.
Vector< Vector< unsigned > > Sparse_assemble_with_arrays_previous_allocation
the number of elements in each row of a compressed matrix in the previous matrix assembly.
virtual void actions_after_parameter_increase(double *const ¶meter_pt)
Empty virtual function; provides hook to perform actions after the increase in the arclength paramete...
friend class PitchForkHandler
void calculate_continuation_derivatives(double *const ¶meter_pt)
A function to calculate the derivatives wrt the arc-length. This version of the function actually doe...
void store_current_dof_values()
Store the current values of the degrees of freedom.
double & dof_derivative(const unsigned &i)
Access function to the derivative of the i-th (local) dof with respect to the arc length,...
bool Problem_has_been_distributed
Has the problem been distributed amongst multiple processors?
void synchronise_all_dofs()
Perform all required synchronisation in solvers.
void get_all_error_estimates(Vector< Vector< double > > &elemental_error)
Return the error estimates computed by (all) refineable (sub)mesh(es) in the elemental_error structur...
bool Empty_actions_before_read_unstructured_meshes_has_been_called
Boolean to indicate that empty actions_before_read_unstructured_meshes() function has been called.
bool Mass_matrix_has_been_computed
Has the mass matrix been computed (and can therefore be reused) Default: false.
unsigned nglobal_data() const
Return the number of global data values.
virtual void actions_before_adapt()
Actions that are to be performed before a mesh adaptation. These might include removing any additiona...
void newton_solve()
Use Newton method to solve the problem.
bool First_jacobian_sign_change
Boolean to indicate whether a sign change has occured in the Jacobian.
void calculate_continuation_derivatives_fd_helper(double *const ¶meter_pt)
A function that performs the guts of the continuation derivative calculation in arc-length continuati...
double Continuation_direction
The direction of the change in parameter that will ensure that a branch is followed in one direction ...
unsigned Parallel_sparse_assemble_previous_allocation
The amount of data allocated during the previous parallel sparse assemble. Used to optimise the next ...
void enable_mass_matrix_reuse()
Enable recycling of the mass matrix in explicit timestepping schemes. Useful for timestepping on fixe...
void steady_newton_solve(unsigned const &max_adapt=0)
Solve a steady problem using adaptive Newton's method, but in the context of an overall unstady probl...
bool Scale_arc_length
Boolean to control whether arc-length should be scaled.
double doubly_adaptive_unsteady_newton_solve_helper(const double &dt, const double &epsilon, const unsigned &max_adapt, const unsigned &suppress_resolve_after_spatial_adapt, const bool &first, const bool &shift=true)
Private helper function that actually performs the unsteady "doubly" adaptive Newton solve....
bool Empty_actions_after_read_unstructured_meshes_has_been_called
Boolean to indicate that empty actions_after_read_unstructured_meshes() function has been called.
virtual void get_residuals(DoubleVector &residuals)
Return the fully-assembled residuals Vector for the problem: Virtual so it can be overloaded in for m...
Vector< GeneralisedElement * > Base_mesh_element_pt
Vector to store the correspondence between a root element and its element number within the global me...
void get_fd_jacobian(DoubleVector &residuals, DenseMatrix< double > &jacobian)
Return the fully-assembled Jacobian and residuals, generated by finite differences.
virtual void sparse_assemble_row_or_column_compressed_with_two_vectors(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
virtual Problem * make_copy()
Make and return a pointer to the copy of the problem. A virtual function that must be filled in by th...
double Maximum_dt
Maximum desired dt.
unsigned long set_timestepper_for_all_data(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data=false)
Set all problem data to have the same timestepper (timestepper_pt) Return the new number of dofs in t...
void activate_pitchfork_tracking(double *const ¶meter_pt, const DoubleVector &symmetry_vector, const bool &block_solve=true)
Turn on pitchfork tracking using the augmented system specified in the PitchForkHandler class....
virtual void dump(std::ofstream &dump_file) const
Dump refinement pattern of all refineable meshes and all generic Problem data to file for restart.
bool Use_finite_differences_for_continuation_derivatives
Boolean to specify which scheme to use to calculate the continuation derivatievs.
bool Arc_length_step_taken
Boolean to indicate whether an arc-length step has been taken.
void set_pinned_values_to_zero()
Set all pinned values to zero. Used to set boundary conditions to be homogeneous in the copy of the p...
bool Default_set_initial_condition_called
Has default set_initial_condition function been called? Default: false.
void get_bifurcation_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction calculated as part of a bifurcation tracking process. If we are not tracking...
double Relaxation_factor
Relaxation fator for Newton method (only a fractional Newton correction is applied if this is less th...
void add_time_stepper_pt(TimeStepper *const &time_stepper_pt)
Add a timestepper to the problem. The function will automatically create or resize the Time object so...
void refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
Refine (one and only!) mesh by splitting the elements identified by their numbers relative to the pro...
LinearAlgebraDistribution *const & dof_distribution_pt() const
Return the pointer to the dof distribution (read-only)
void prune_halo_elements_and_nodes(DocInfo &doc_info, const bool &report_stats)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement,...
virtual void get_inverse_mass_matrix_times_residuals(DoubleVector &Mres)
Return the residual vector multiplied by the inverse mass matrix Virtual so that it can be overloaded...
void check_halo_schemes()
Check the halo/haloed node/element schemes.
double Parameter_current
Storage for the present value of the global parameter.
double *& dof_pt(const unsigned &i)
Pointer to i-th dof in the problem.
void assign_eigenvector_to_dofs(DoubleVector &eigenvector)
Assign the eigenvector passed to the function to the dofs in the problem so that it can be output by ...
@ Uniform_matrix_distribution
@ Default_matrix_distribution
@ Problem_matrix_distribution
void send_refinement_info_helper(Vector< unsigned > &old_domain_for_base_element, Vector< unsigned > &new_domain_for_base_element, const unsigned &max_refinement_level_overall, std::map< unsigned, Vector< unsigned > > &refinement_info_for_root_local, Vector< Vector< Vector< unsigned > > > &refinement_info_for_root_elements)
Send refinement information between processors.
unsigned setup_element_count_per_dof()
Function that populates the Element_counter_per_dof vector with the number of elements that contribut...
OomphCommunicator * Communicator_pt
The communicator for this problem.
double Newton_solver_tolerance
The Tolerance below which the Newton Method is deemed to have converged.
void set_consistent_pinned_values_for_continuation()
Private helper function that is used to set the appropriate pinned values for continuation.
void activate_bifurcation_tracking(double *const ¶meter_pt, const DoubleVector &eigenvector, const bool &block_solve=true)
Activate generic bifurcation tracking for a single (real) eigenvalue where the initial guess for the ...
LinearSolver *& mass_matrix_solver_for_explicit_timestepper_pt()
Return a pointer to the linear solver object used for explicit time stepping.
bool Discontinuous_element_formulation
Is the problem a discontinuous one, i.e. can the elemental contributions be treated independently....
bool Doc_imbalance_in_parallel_assembly
Boolean to switch on assessment of load imbalance in parallel assembly of distributed problem.
long synchronise_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Classify any non-classified nodes into halo/haloed and synchronise equation numbers....
void create_new_linear_algebra_distribution(LinearAlgebraDistribution *&dist_pt)
Get new linear algebra distribution (you're in charge of deleting it!)
void p_refine_uniformly_aux(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info, const bool &prune)
Helper function to do compund p-refinement of (all) p-refineable (sub)mesh(es) uniformly as many time...
void calculate_continuation_derivatives_helper(const DoubleVector &z)
A function that performs the guts of the continuation derivative calculation in arc length continuati...
Vector< Problem * > Copy_of_problem_pt
Vector of pointers to copies of the problem used in adaptive bifurcation tracking problems (ALH: TEMP...
Vector< unsigned > distribute(const Vector< unsigned > &element_partition, DocInfo &doc_info, const bool &report_stats=false)
Distribute the problem and doc, using the specified partition; returns a vector which details the par...
Mesh * Mesh_pt
The mesh pointer.
double Parameter_derivative
Storage for the derivative of the global parameter wrt arc-length.
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
void add_eigenvector_to_dofs(const double &epsilon, const DoubleVector &eigenvector)
Add the eigenvector passed to the function scaled by the constat epsilon to the dofs in the problem s...
Vector< Data * > Global_data_pt
Vector of global data: "Nobody" (i.e. none of the elements etc.) is "in charge" of this Data so it wo...
void recompute_load_balanced_assembly()
Helper function to re-assign the first and last elements to be assembled by each processor during par...
Vector< double * > Dof_pt
Vector of pointers to dofs.
void p_adapt()
p-adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on their own error...
double DTSF_max_increase
Maximum possible increase of dt between time-steps in adaptive schemes.
bool Must_recompute_load_balance_for_assembly
Boolean indicating that the division of elements over processors during the assembly process must be ...
virtual void shift_time_values()
Shift all values along to prepare for next timestep.
double Target_error_safety_factor
Safety factor to ensure we are aiming for a target error, TARGET, that is below our tolerance: TARGET...
void set_default_first_and_last_element_for_assembly()
Set default first and last elements for parallel assembly of non-distributed problem.
bool Keep_temporal_error_below_tolerance
Boolean to decide if a timestep is to be rejected if the error estimate post-solve (computed by globa...
bool Shut_up_in_newton_solve
Boolean to indicate if all output is suppressed in Problem::newton_solve(). Defaults to false.
void set_dofs(const DoubleVector &dofs)
Set the values of the dofs.
void setup_base_mesh_info_after_pruning()
Helper function to re-setup the Base_mesh enumeration (used during load balancing) after pruning.
Vector< Mesh * > Sub_mesh_pt
Vector of pointers to submeshes.
ExplicitTimeStepper *& explicit_time_stepper_pt()
Return a pointer to the explicit timestepper.
Vector< double > Max_res
Maximum residuals at start and after each newton iteration.
EigenSolver * Eigen_solver_pt
Pointer to the eigen solver for the problem.
bool Bisect_to_find_bifurcation
Boolean to control wheter bisection is used to located bifurcation.
void explicit_timestep(const double &dt, const bool &shift_values=true)
Take an explicit timestep of size dt and optionally shift any stored values of the time history.
void delete_all_external_storage()
Wrapper function to delete external storage for each submesh of the problem.
friend class BlockPitchForkLinearSolver
double DTSF_min_decrease
Minimum allowed decrease of dt between time-steps in adaptive schemes. Lower scaling values will reje...
virtual void symmetrise_eigenfunction_for_adaptive_pitchfork_tracking()
Virtual function that is used to symmetrise the problem so that the current solution exactly satisfie...
double Minimum_dt_but_still_proceed
If Minimum_dt_but_still_proceed positive, then dt will not be reduced below this value during adaptiv...
double adaptive_unsteady_newton_solve(const double &dt_desired, const double &epsilon)
Attempt to advance timestep by dt_desired. If the solution fails the timestep will be halved until co...
unsigned Desired_newton_iterations_ds
The desired number of Newton Steps to reach convergence at each step along the arc.
void rebuild_global_mesh()
If one of the submeshes has changed (e.g. by mesh adaptation) we need to update the global mesh....
void solve_adjoint_eigenproblem(const unsigned &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector_real, Vector< DoubleVector > &eigenvector_imag, const bool &steady=true)
Solve an adjoint eigenvalue problem using the same procedure as solve_eigenproblem....
void activate_hopf_tracking(double *const ¶meter_pt, const bool &block_solve=true)
Turn on Hopf bifurcation tracking using the augmented system specified in the HopfHandler class....
Vector< double * > Halo_dof_pt
Storage for the halo degrees of freedom (only required) when accessing via the global equation number...
void deactivate_bifurcation_tracking()
Deactivate all bifuraction tracking, by reseting the assembly handler to the default.
void globally_convergent_line_search(const Vector< double > &x_old, const double &half_residual_squared_old, DoubleVector &gradient, DoubleVector &newton_dir, double &half_residual_squared, const double &stpmax)
Line search helper for globally convergent Newton method.
void get_hessian_vector_products(DoubleVectorWithHaloEntries const &Y, Vector< DoubleVectorWithHaloEntries > const &C, Vector< DoubleVectorWithHaloEntries > &product)
Return the product of the global hessian (derivative of Jacobian matrix with respect to all variables...
virtual double global_temporal_error_norm()
Function to calculate a global error norm, used in adaptive timestepping to control the change in tim...
@ Perform_assembly_using_two_arrays
@ Perform_assembly_using_maps
@ Perform_assembly_using_two_vectors
@ Perform_assembly_using_vectors_of_pairs
@ Perform_assembly_using_lists
void refine_distributed_base_mesh(Vector< Vector< Vector< unsigned > > > &to_be_refined_on_each_root, const unsigned &max_level_overall)
Load balance helper routine: refine each new base (sub)mesh based upon the elements to be refined wit...
int Sign_of_jacobian
Storage for the sign of the global Jacobian.
std::map< GeneralisedElement *, unsigned > Base_mesh_element_number_plus_one
Map which stores the correspondence between a root element and its element number (plus one) within t...
double Max_residuals
Maximum desired residual: if the maximum residual exceeds this value, the program will exit.
unsigned nsub_mesh() const
Return number of submeshes.
double & time()
Return the current value of continuous time.
unsigned self_test()
Self-test: Check meshes and global data. Return 0 for OK.
unsigned ntime_stepper() const
Return the number of time steppers.
void activate_fold_tracking(double *const ¶meter_pt, const bool &block_solve=true)
Turn on fold tracking using the augmented system specified in the FoldHandler class....
bool Use_globally_convergent_newton_method
Use the globally convergent newton method.
LinearSolver * Default_linear_solver_pt
Pointer to the default linear solver.
double Minimum_ds
Minimum desired value of arc-length.
void get_data_to_be_sent_during_load_balancing(const Vector< unsigned > &element_domain_on_this_proc, Vector< int > &send_n, Vector< double > &send_data, Vector< int > &send_displacement, Vector< unsigned > &old_domain_for_base_element, Vector< unsigned > &new_domain_for_base_element, unsigned &max_refinement_level_overall)
Load balance helper routine: Get data to be sent to other processors during load balancing and other ...
void restore_dof_values()
Restore the stored values of the degrees of freedom.
double Numerical_zero_for_sparse_assembly
A tolerance used to determine whether the entry in a sparse matrix is zero. If it is then storage nee...
double FD_step_used_in_get_hessian_vector_products
virtual void partition_global_mesh(Mesh *&global_mesh_pt, DocInfo &doc_info, Vector< unsigned > &element_domain, const bool &report_stats=false)
Partition the global mesh, return vector specifying the processor number for each element....
unsigned Sparse_assemble_with_arrays_initial_allocation
the number of elements to initially allocate for a matrix row within the sparse_assembly_with_two_arr...
void bifurcation_adapt_doc_errors(const unsigned &bifurcation_type)
A function that is used to document the errors used in the adaptive solution of bifurcation problems.
double Ds_current
Storage for the current step value.
virtual void set_initial_condition()
Set initial condition (incl previous timesteps). We need to establish this interface because I....
LinearSolver * Mass_matrix_solver_for_explicit_timestepper_pt
Pointer to the linear solver used for explicit time steps (this is likely to be different to the line...
void get_all_halo_data(std::map< unsigned, double * > &map_of_halo_data)
Get pointers to all possible halo data indexed by global equation number in a map.
void load_balance()
Balance the load of a (possibly non-uniformly refined) problem that has already been distributed,...
double Minimum_dt
Minimum desired dt: if dt falls below this value, exit.
double arc_length_step_solve(double *const ¶meter_pt, const double &ds, const unsigned &max_adapt=0)
Solve a steady problem using arc-length continuation, when the parameter that becomes a variable corr...
virtual void actions_after_adapt()
Actions that are to be performed after a mesh adaptation.
void p_refine_uniformly()
p-refine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem
bool Use_continuation_timestepper
Boolean to control original or new storage of dof stuff.
void calculate_continuation_derivatives_fd(double *const ¶meter_pt)
A function to calculate the derivatives with respect to the arc-length required for continuation by f...
LinearAlgebraDistribution * Dof_distribution_pt
The distribution of the DOFs in this problem. This object is created in the Problem constructor and s...
void refine_uniformly()
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem.
static ContinuationStorageScheme Continuation_time_stepper
Storage for the single static continuation timestorage object.
bool Problem_is_nonlinear
Boolean flag indicating if we're dealing with a linear or nonlinear Problem – if set to false the New...
virtual void sparse_assemble_row_or_column_compressed_with_maps(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
AssemblyHandler * Default_assembly_handler_pt
Pointer to the default assembly handler.
void get_my_eqns(AssemblyHandler *const &assembly_handler_pt, const unsigned &el_lo, const unsigned &el_hi, Vector< unsigned > &my_eqns)
Helper method that returns the (unique) global equations to which the elements in the range el_lo to ...
virtual void read(std::ifstream &restart_file, bool &unsteady_restart)
Read refinement pattern of all refineable meshes and refine them accordingly, then read all Data and ...
virtual ~Problem()
Virtual destructor to clean up memory.
Mesh *& mesh_pt()
Return a pointer to the global mesh.
void get_dofs(DoubleVector &dofs) const
Return the vector of dofs, i.e. a vector containing the current values of all unknowns.
Time * Time_pt
Pointer to global time for the problem.
DoubleVectorWithHaloEntries Element_count_per_dof
Counter that records how many elements contribute to each dof. Used to determine the (discrete) arc-l...
virtual void actions_before_newton_convergence_check()
Any actions that are to be performed before the residual is checked in the Newton method,...
Time *& time_pt()
Return a pointer to the global time object.
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
unsigned newton_solve_continuation(double *const ¶meter_pt)
Perform a basic arc-length continuation step using Newton's method. Returns number of Newton steps ta...
double Max_permitted_error_for_halo_check
Threshold for error throwing in Problem::check_halo_schemes()
unsigned Sparse_assemble_with_arrays_allocation_increment
the number of elements to add to a matrix row when the initial allocation is exceeded within the spar...
void reset_assembly_handler_to_default()
Reset the system to the standard non-augemented state.
void solve_eigenproblem(const unsigned &n_eval, Vector< std::complex< double > > &alpha, Vector< double > &beta, Vector< DoubleVector > &eigenvector_real, Vector< DoubleVector > &eigenvector_imag, const bool &steady=true)
Solve an eigenproblem as assembled by the Problem's constituent elements. Calculate (at least) n_eval...
virtual void actions_after_newton_step()
Any actions that are to be performed after each individual Newton step. Most likely to be used for di...
bool Pause_at_end_of_sparse_assembly
Protected boolean flag to halt program execution during sparse assemble process to assess peak memory...
void remove_null_pointers_from_external_halo_node_storage()
Consolidate external halo node storage by removing nulled out pointers in external halo and haloed sc...
Vector< double > * Saved_dof_pt
Pointer to vector for backup of dofs.
void unsteady_newton_solve(const double &dt)
Advance time by dt and solve by Newton's method. This version always shifts time values.
AssemblyHandler * Assembly_handler_pt
virtual void actions_before_read_unstructured_meshes()
Actions that are to be performed before reading in restart data for problems involving unstructured b...
Vector< TimeStepper * > Time_stepper_pt
The Vector of time steppers (there could be many different ones in multiphysics problems)
void doc_errors()
Get max and min error for all elements in submeshes.
bool Mass_matrix_reuse_is_enabled
Is re-use of the mass matrix in explicit timestepping enabled Default:false.
void get_derivative_wrt_global_parameter(double *const ¶meter_pt, DoubleVector &result)
Get the derivative of the entire residuals vector wrt a global parameter, used in continuation proble...
bool distributed() const
If we have MPI return the "problem has been distributed" flag, otherwise it can't be distributed so r...
bool are_hessian_products_calculated_analytically()
Function to determine whether the hessian products are calculated analytically.
bool does_pointer_correspond_to_problem_data(double *const ¶meter_pt)
Return a boolean flag to indicate whether the pointer parameter_pt refers to values stored in a Data ...
ExplicitTimeStepper * Explicit_time_stepper_pt
Pointer to a single explicit timestepper.
double arc_length_step_solve_helper(double *const ¶meter_pt, const double &ds, const unsigned &max_adapt)
Private helper function that actually contains the guts of the arc-length stepping,...
bool Use_predictor_values_as_initial_guess
Use values from the time stepper predictor as an initial guess.
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
Base class for refineable meshes. Provides standardised interfaces for the following standard mesh ad...
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
SuperLU Project Solver class. This is a combined wrapper for both SuperLU and SuperLU Dist....
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
virtual unsigned ndt() const =0
Number of timestep increments that are required by the scheme.
virtual void shift_time_values(Data *const &data_pt)=0
This function advances the Data's time history so that we can move on to the next timestep.
virtual void set_weights()=0
Function to set the weights for present timestep (don't need to pass present timestep or previous tim...
ExplicitTimeStepper * explicit_predictor_pt()
Get the pointer to the explicit timestepper to use as a predictor in adaptivity if Predict_by_explici...
virtual void calculate_predicted_values(Data *const &data_pt)
Do the predictor step for data stored in a Data object (currently empty – overwrite for specific sche...
virtual void set_predictor_weights()
Set the weights for the predictor previous timestep (currently empty – overwrite for specific scheme)
virtual void actions_before_timestep(Problem *problem_pt)
Interface for any actions that need to be performed before a time step.
virtual void actions_after_timestep(Problem *problem_pt)
Interface for any actions that need to be performed after a time step.
virtual void assign_initial_values_impulsive(Data *const &data_pt)=0
Initialise the time-history for the Data values corresponding to an impulsive start.
void make_steady()
Function to make the time stepper temporarily steady. This is trivially achieved by setting all the w...
virtual void undo_make_steady()
Reset the is_steady status of a specific TimeStepper to its default and re-assign the weights.
void update_predicted_time(const double &new_time)
Set the time that the current predictions apply for, only needed for paranoid checks when doing Predi...
Time *const & time_pt() const
Access function for the pointer to time (const version)
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
virtual void set_error_weights()
Set the weights for the error computation, (currently empty – overwrite for specific scheme)
Class to keep track of discrete/continous time. It is essential to have a single Time object when usi...
double & time()
Return the current value of the continuous time.
void shift_dt()
Update all stored values of dt by shifting each value along the array. This function must be called b...
void initialise_dt(const double &dt_)
Set all timesteps to the same value, dt.
double & dt(const unsigned &t=0)
Return the value of the t-th stored timestep (t=0: present; t>0: previous).
unsigned ndt() const
Return the number of timesteps stored.
void resize(const unsigned &n_dt)
Resize the vector holding the number of previous timesteps and initialise the new values to zero.
Base class for tree-based refineable meshes.
TreeRoot is a Tree that forms the root of a (recursive) tree. The "root node" is special as it holds ...
Base class for triangle meshes (meshes made of 2D triangle elements). Note: we choose to template Tri...
A slight extension to the standard template vector class so that we can include "graceful" array rang...
bool Doc_comprehensive_timings
Global boolean to switch on comprehensive timing – can probably be declared const false when developm...
void partition_distributed_mesh(Problem *problem_pt, const unsigned &objective, Vector< unsigned > &element_domain_on_this_proc, const bool &bypass_metis=false)
Use METIS to assign each element in an already-distributed mesh to a domain. On return,...
void partition_mesh(Problem *problem_pt, const unsigned &ndomain, const unsigned &objective, Vector< unsigned > &element_domain)
Use METIS to assign each element to a domain. On return, element_domain[ielem] contains the number of...
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types).
void clean_up_memory()
Clean up function that deletes anything dynamically allocated in this namespace.
void setup()
Setup terminate helper.
double timer()
returns the time in seconds after some point in past
std::string convert_secs_to_formatted_string(const double &time_in_sec)
Returns a nicely formatted string from an input time in seconds; the format depends on the size of ti...
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
void pause(std::string message)
Pause and display message.
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...