72 Explicit_time_stepper_pt(0),
74 Default_set_initial_condition_called(
false),
75 Use_globally_convergent_newton_method(
false),
76 Empty_actions_before_read_unstructured_meshes_has_been_called(
false),
77 Empty_actions_after_read_unstructured_meshes_has_been_called(
false),
78 Store_local_dof_pt_in_elements(
false),
79 Calculate_hessian_products_analytic(
false),
81 Doc_imbalance_in_parallel_assembly(
false),
82 Use_default_partition_in_load_balance(
false),
83 Must_recompute_load_balance_for_assembly(
true),
86 Relaxation_factor(1.0),
87 Newton_solver_tolerance(1.0e-8),
88 Max_newton_iterations(10),
89 Nnewton_iter_taken(0),
91 Time_adaptive_newton_crash_on_solve_fail(
false),
92 Jacobian_reuse_is_enabled(
false),
93 Jacobian_has_been_computed(
false),
94 Problem_is_nonlinear(
true),
95 Pause_at_end_of_sparse_assembly(
false),
96 Doc_time_in_distribute(
false),
97 Sparse_assembly_method(Perform_assembly_using_vectors_of_pairs),
98 Sparse_assemble_with_arrays_initial_allocation(400),
99 Sparse_assemble_with_arrays_allocation_increment(150),
100 Numerical_zero_for_sparse_assembly(0.0),
101 FD_step_used_in_get_hessian_vector_products(1.0e-8),
102 Mass_matrix_reuse_is_enabled(
false),
103 Mass_matrix_has_been_computed(
false),
104 Discontinuous_element_formulation(
false),
107 DTSF_max_increase(4.0),
108 DTSF_min_decrease(0.8),
109 Target_error_safety_factor(1.0),
110 Minimum_dt_but_still_proceed(-1.0),
111 Scale_arc_length(
true),
112 Desired_proportion_of_arc_length(0.5),
115 Continuation_direction(1.0),
116 Parameter_derivative(1.0),
117 Parameter_current(0.0),
118 Use_continuation_timestepper(
false),
119 Dof_derivative_offset(1),
120 Dof_current_offset(2),
122 Desired_newton_iterations_ds(5),
124 Bifurcation_detection(
false),
125 Bisect_to_find_bifurcation(
false),
126 First_jacobian_sign_change(
false),
127 Arc_length_step_taken(
false),
128 Use_finite_differences_for_continuation_derivatives(
false),
130 Dist_problem_matrix_distribution(Uniform_matrix_distribution),
131 Parallel_sparse_assemble_previous_allocation(0),
132 Problem_has_been_distributed(
false),
133 Bypass_increase_in_dof_check_during_pruning(
false),
134 Max_permitted_error_for_halo_check(1.0e-14),
136 Shut_up_in_newton_solve(
false),
137 Always_take_one_newton_step(
false),
138 Timestep_reduction_factor_after_nonconvergence(0.5),
139 Keep_temporal_error_below_tolerance(
true)
208 for (
unsigned c = 0; c <
n_copies; c++)
257 for (
unsigned n = 0;
n <
n_var;
n++)
268 unsigned Nelement = 0;
304 const unsigned nrow = this->
ndof();
344 for (
unsigned p = 0;
p <
nproc;
p++)
368 error_stream <<
"Never get here. Dist_problem_matrix_distribution = "
442 for (
unsigned i = 0;
i <
n;
i++)
453 warn_message <<
"WARNING: All entries in specified partitioning vector \n"
454 <<
" are zero -- will ignore this and use METIS\n"
455 <<
" to perform the partitioning\n";
510 warn_message <<
"WARNING: You've tried to distribute a problem over\n"
511 <<
"only one processor: this would make METIS crash.\n"
512 <<
"Ignoring your request for distribution.\n";
514 "Problem::distribute()",
522 error_stream <<
"You have tried to distribute a problem\n"
523 <<
"but there are less elements than processors.\n"
524 <<
"Please re-run with more elements!\n"
525 <<
"Please also ensure that actions_before_distribute().\n"
526 <<
"and actions_after_distribute() are correctly set up.\n"
581 error_stream <<
"You have tried to distribute a problem\n"
582 <<
"but at least one of your meshes is no longer\n"
583 <<
"uniformly refined. In order to preserve the Tree\n"
584 <<
"and TreeForest structure, Problem::distribute() can\n"
585 <<
"only be called while meshes are uniformly refined.\n"
596 error_stream <<
"You have tried to distribute a problem\n"
597 <<
"and there is some global data.\n"
598 <<
"This is not likely to work...\n"
644 oomph_info <<
"INFO: using METIS to partition elements" << std::endl;
650 oomph_info <<
"INFO: using pre-set partition of elements"
663 oomph_info <<
"Time for partitioning of global mesh: "
687 oomph_info <<
"Time for actions before distribute: "
756 for (
unsigned e = 0;
e <
n_ele;
e++)
806 for (
unsigned e = 0;
e <
n_ele;
e++)
823 <<
") is not the sameas the number of\nadded elements ("
824 <<
counter <<
") to the Base_mesh_element_pt data "
825 <<
"structure!!!\n\n";
827 "Problem::distribute()",
842 oomph_info <<
"INFO: We're over-ruling the \"keep as halo element\"\n"
843 <<
" status because the preset partitioning\n"
844 <<
" didn't place ANY elements on this processor,\n"
845 <<
" probably because of a restart on a larger \n"
846 <<
" number of processors\n";
869 <<
"--------------------" << std::endl;
874 submesh_element_domain[
i_mesh],
886 for (
unsigned e = 0;
e <
n_del;
e++)
923 oomph_info <<
"Time for re-assigning eqn numbers (in distribute): "
933 <<
"Number of dofs in distribute() has changed "
935 <<
"Check that you've implemented any necessary "
936 "actions_before/after\n"
937 <<
"distribute functions, e.g. to pin redundant pressure dofs"
984 << doc_info.
number() <<
".dat";
1015 for (
unsigned e = 0;
e <
nelem;
e++)
1024 for (
unsigned e = 0;
e <
nelem;
e++)
1047 for (
int d = 0; d <
n_proc; d++)
1069 error_stream <<
"No process has more than 1 element, and\n"
1070 <<
"at least one process has no elements!\n"
1071 <<
"Suggest rerunning with more refinement.\n"
1080 for (
unsigned e = 0;
e <
nelem;
e++)
1093 oomph_info <<
"INFO: Switched element domain at position " <<
e
1096 <<
" to process " << d << std::endl
1097 <<
"which was given no elements by METIS partition"
1110 for (
unsigned e = 0;
e <
nelem;
e++)
1116 for (
unsigned e = 0;
e <
nelem;
e++)
1127 <<
" elements from this partition" << std::endl
1150 <<
"WARNING: Problem::prune_halo_elements_and_nodes() was called on a "
1151 <<
"non-distributed Problem!" << std::endl;
1152 oomph_info <<
"Ignoring your request..." << std::endl;
1160 <<
"WARNING: You've tried to prune halo layers on a problem\n"
1161 <<
"with only one processor: this is unnecessary.\n"
1162 <<
"Ignoring your request." << std::endl
1184 oomph_info <<
"Time for actions_before_distribute() in "
1185 <<
"Problem::prune_halo_elements_and_nodes(): "
1192 std::map<GeneralisedElement*, unsigned>
1195 for (
unsigned e = 0;
e <
nel;
e++)
1220 ref_el_pt->tree_pt()->stick_all_tree_nodes_into_vector(tree_pt);
1221 unsigned ntree = tree_pt.
size();
1222 for (
unsigned t = 0;
t < ntree;
t++)
1235 oomph_info <<
"Time for setup old root elements in "
1236 <<
"Problem::prune_halo_elements_and_nodes(): "
1271 oomph_info <<
"Total time for all mesh-level prunes in "
1272 <<
"Problem::prune_halo_elements_and_nodes(): "
1364 ref_el_pt->tree_pt()->root_pt()->object_pt();
1401 oomph_info <<
"Time for setup of new base elements in "
1402 <<
"Problem::prune_halo_elements_and_nodes(): "
1423 <<
"Problem::prune_halo_elements_and_nodes(): "
1448 <<
"Number of new root elements spawned from old root " <<
e
1508 oomph_info <<
"Time for finishing off base mesh info "
1509 <<
"Problem::prune_halo_elements_and_nodes(): "
1522 oomph_info <<
"Time for actions_after_distribute() "
1523 <<
"Problem::prune_halo_elements_and_nodes(): "
1540 oomph_info <<
"Time for assign_eqn_numbers() "
1541 <<
"Problem::prune_halo_elements_and_nodes(): "
1554 <<
"Number of dofs in prune_halo_elements_and_nodes() has "
1557 <<
"Check that you've implemented any necessary "
1558 "actions_before/after"
1559 <<
"\nadapt/distribute functions, e.g. to pin redundant pressure"
1586 std::string error_message =
"Problem::build_global_mesh() called,\n";
1587 error_message +=
" but a global mesh has already been built:\n";
1588 error_message +=
"Problem::Mesh_pt is not zero!\n";
1596 std::string error_message =
"Problem::build_global_mesh() called,\n";
1597 error_message +=
" but there are no submeshes:\n";
1598 error_message +=
"Problem::Sub_mesh_pt has no entries\n";
1645 oomph_info <<
"Created Time with " << ndt <<
" timesteps" << std::endl;
1654 oomph_info <<
"Resized Time to include " << ndt <<
" timesteps"
1660 oomph_info <<
"Time object already has storage for " << ndt
1661 <<
" timesteps" << std::endl;
1684 oomph_info <<
"Created Time with storage for no previous timestep"
1689 oomph_info <<
"Time object already exists " << std::endl;
1750 oomph_info <<
"Problem is not distributed. Parallel assembly of "
1751 <<
"Jacobian uses default partitioning: " << std::endl;
1756 oomph_info <<
"Proc " <<
p <<
" assembles from element "
1762 oomph_info <<
"Proc " <<
p <<
" assembles no elements\n";
1790 oomph_info <<
"Not re-computing distribution of elemental assembly\n"
1791 <<
"because there are fewer elements than processors\n";
1815 for (
unsigned e = 0;
e <
nel;
e++)
1848 <<
"Re-assigning distribution of element assembly over processors:"
1919 error_stream <<
"Error: First/last element of proc " <<
p <<
": "
1927 error_stream <<
"Error: First element of proc " <<
p + 1 <<
": "
2009 oomph_info <<
"Processor " << 0 <<
" assembles Jacobians"
2028 oomph_info <<
"Processor " <<
p <<
" assembles Jacobians"
2077 const bool& assign_local_eqn_numbers)
2084 error_stream <<
"Global mesh does not exist, so equation numbers cannot "
2089 error_stream <<
"There aren't even any sub-meshes in the Problem.\n"
2090 <<
"You can set the global mesh directly by using\n"
2091 <<
"Problem::mesh_pt() = my_mesh_pt;\n"
2092 <<
"OR you can use Problem::add_sub_mesh(mesh_pt); "
2093 <<
"to add a sub mesh.\n";
2099 error_stream <<
"You need to call Problem::build_global_mesh() to create "
2101 <<
"from the sub-meshes.\n\n";
2152 for (
unsigned e = 0;
e <
nel;
e++)
2181 <<
"Time for complete setup of dependencies in assign_eqn_numbers: "
2226 n_dof = spine_mesh_pt->assign_global_spine_eqn_numbers(
Dof_pt);
2238 n_dof = spine_mesh_pt->assign_global_spine_eqn_numbers(
Dof_pt);
2247 <<
"Time for assign_global_eqn_numbers in assign_eqn_numbers: "
2280 oomph_info <<
"Time for Problem::synchronise_eqn_numbers in "
2324 std::stringstream
tmp;
2325 tmp <<
"Time for calls to Problem::remove_duplicate_data in "
2332 tmp <<
" removed some/any data.\n";
2354 std::stringstream
tmp;
2356 <<
"Time for MPI_Allreduce after Problem::remove_duplicate_data in "
2357 <<
"Problem::assign_eqn_numbers: " <<
t_end -
t_start << std::endl;
2379 <<
"Problem::remove_null_pointers_from_external_halo_node_"
2411 if (assign_local_eqn_numbers)
2430 oomph_info <<
"Total time for all Mesh::assign_local_eqn_numbers in "
2454 <<
"Global mesh does not exist, so equation numbers cannot be found.\n";
2458 error_stream <<
"There aren't even any sub-meshes in the Problem.\n"
2459 <<
"You can set the global mesh directly by using\n"
2460 <<
"Problem::mesh_pt() = my_mesh_pt;\n"
2461 <<
"OR you can use Problem::add_sub_mesh(mesh_pt); "
2462 <<
"to add a sub mesh.\n";
2468 error_stream <<
"You need to call Problem::build_global_mesh() to create "
2470 <<
"from the sub-meshes.\n\n";
2478 <<
"Although this program will describe the degrees of freedom in the \n"
2479 <<
"problem, it will do so using the typedef for the elements. This is \n"
2480 <<
"not neccesarily human readable, but there is a solution.\n"
2481 <<
"Pipe your program's output through c++filt, with the argument -t.\n"
2482 <<
"e.g. \"./two_d_multi_poisson | c++filt -t > ReadableOutput.txt\".\n "
2483 <<
"(Disregarding the quotes)\n\n\n";
2485 out <<
"Classifying Global Equation Numbers" << std::endl;
2486 out << std::string(80,
'-') << std::endl;
2504 std::string
in(
" in Problem's Only Mesh.");
2514 std::string
in(
" in Problem's Only SpineMesh.");
2515 spine_mesh_pt->describe_spine_dofs(
out,
in);
2530 spine_mesh_pt->describe_spine_dofs(
out,
in);
2536 out << std::string(80,
'\\') << std::endl;
2537 out << std::string(80,
'\\') << std::endl;
2538 out << std::string(80,
'\\') << std::endl;
2539 out <<
"Classifying global eqn numbers in terms of elements." << std::endl;
2540 out << std::string(80,
'-') << std::endl;
2541 out <<
"Eqns | Source" << std::endl;
2542 out << std::string(80,
'-') << std::endl;
2546 std::string
in(
" in Problem's Only Mesh.");
2575 for (
unsigned long l = 0;
l <
n_dof;
l++)
2587 throw OomphLibError(
"Not designed for distributed problems",
2605 if (eqn_number >= 0)
2622 if (eqn_number >= 0)
2638 if (eqn_number >= 0)
2670 actually_removed_some_data =
false;
2761 ->ncont_interpolated_values();
2930 unsigned nb = (*boundaries_pt).
size();
2932 for (std::set<unsigned>::iterator
it =
2933 (*boundaries_pt).begin();
2934 it != (*boundaries_pt).end();
2939 for (
unsigned i = 0;
i <
nb;
i++)
2968 ->ncont_interpolated_values();
2995 <<
"Number of master nodes for node to be replaced, "
2998 <<
" for i_cont=" <<
i_cont << std::endl;
3001 <<
"Nodal coordinates of replacement node:";
3003 for (
unsigned i = 0;
i < ndim;
i++)
3010 <<
" master nodes are: \n";
3015 ->master_node_pt(
k);
3017 for (
unsigned i = 0;
i < ndim;
i++)
3027 <<
"Nodal coordinates of node to be replaced:";
3029 for (
unsigned i = 0;
i < ndim;
i++)
3037 <<
" master nodes are: \n";
3044 for (
unsigned i = 0;
i < ndim;
i++)
3079 ->ncont_interpolated_values();
3157 for (std::set<unsigned>::iterator
it =
3158 (*boundaries_pt).begin();
3159 it != (*boundaries_pt).end();
3186 <<
"About to re-set master for i_cont= " <<
i_cont
3187 <<
" for external node (with proc " <<
iproc
3190 for (
unsigned jj = 0;
jj <
n;
jj++)
3195 <<
" which is not hanging --> About to die!"
3196 <<
"Outputting offending element into oomph-info "
3333 for (
unsigned j = 0;
j <
nnod;
j++)
3470 for (
unsigned j = 0;
j <
nnod;
j++)
3500 const unsigned long n_dof = this->
ndof();
3502 if (n_dof !=
dofs.nrow())
3505 error_stream <<
"Number of degrees of freedom in vector argument "
3506 <<
dofs.nrow() <<
"\n"
3507 <<
"does not equal number of degrees of freedom in problem "
3513 for (
unsigned long l = 0;
l <
n_dof;
l++)
3525 throw OomphLibError(
"Not designed for distributed problems",
3540 if (eqn_number >= 0)
3557 if (eqn_number >= 0)
3573 if (eqn_number >= 0)
3589 throw OomphLibError(
"Not implemented for distributed problems!",
3607 if (eqn_number >= 0)
3623 if (eqn_number >= 0)
3639 if (eqn_number >= 0)
3654 const unsigned long n_dof = this->
ndof();
3655 for (
unsigned long l = 0;
l <
n_dof;
l++)
3676 error_stream <<
"The function get_inverse_mass_matrix_times_residuals() "
3678 <<
"used with the default assembly handler\n\n";
3725 oomph_info <<
"Not recomputing Mass Matrix " << std::endl;
3744 oomph_info <<
"Enabling resolve in explicit timestep" << std::endl;
3823 error_stream <<
"The distribution of the residuals vector does not "
3824 "have the correct\n"
3825 <<
"number of global rows\n";
3899 dist_pt, column_index, row_start, value, nnz,
res);
3934 if (
residuals.distribution_pt()->distributed())
3938 <<
"If the DoubleVector residuals is setup then it must not "
3939 <<
"be distributed.";
3947 <<
"If the DoubleVector residuals is setup then it must have"
3948 <<
" the correct number of rows";
3953 *
residuals.distribution_pt()->communicator_pt()))
3957 <<
"If the DoubleVector residuals is setup then it must have"
3958 <<
" the same communicator as the problem.";
4063 error_stream <<
"The distribution of the residuals must "
4064 <<
"be the same as the distribution of the jacobian.";
4072 <<
"The distribution of the jacobian and residuals does not"
4073 <<
"have the correct number of global rows.";
4081 error_stream <<
"The distribution of the jacobian and residuals must "
4082 <<
"both be setup or both not setup";
4119 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4129 dist_pt, column_index, row_start, value, nnz,
res);
4132 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4144 dist_pt->nrow(), nnz[0], value[0], column_index[0], row_start[0]);
4189 if (
residuals.distribution_pt()->distributed())
4193 <<
"If the DoubleVector residuals is setup then it must not "
4194 <<
"be distributed.";
4202 <<
"If the DoubleVector residuals is setup then it must have"
4203 <<
" the correct number of rows";
4208 *
residuals.distribution_pt()->communicator_pt()))
4212 <<
"If the DoubleVector residuals is setup then it must have"
4213 <<
" the same communicator as the problem.";
4256 value[0], row_index[0], column_start[0], nnz[0],
n_dof,
n_dof);
4264 error_stream <<
"Cannot assemble a CCDoubleMatrix Jacobian on more "
4265 <<
"than one processor.";
4336 for (
unsigned i = 0;
i <
n_dim;
i++)
4404 for (
unsigned i = 0;
i <
n_dim;
i++)
4536 <<
"Error: Incorrect value for Problem::Sparse_assembly_method"
4538 <<
"It should be one of the enumeration Problem::Assembly_method"
4574 unsigned long el_lo = 0;
4617 <<
"row_or_column_start.size() "
4619 <<
"column_or_row_index.size() "
4629 <<
"Error in Problem::sparse_assemble_row_or_column_compressed "
4631 <<
"value.size() " << value.
size() <<
" does not equal "
4672 for (
unsigned i = 0;
i <
ndof;
i++)
4742 for (
unsigned i = 0;
i <
nvar;
i++)
4755 for (
unsigned j = 0;
j <
nvar;
j++)
4850 value[
m] =
new double[nnz[
m]];
4866 for (std::map<unsigned, double>::iterator
it =
4886 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
4887 pause(
"Check memory usage now.");
4920 unsigned long el_lo = 0;
4963 <<
"row_or_column_start.size() "
4965 <<
"column_or_row_index.size() "
4975 <<
"Error in Problem::sparse_assemble_row_or_column_compressed "
4977 <<
"value.size() " << value.
size() <<
" does not equal "
5016 for (
unsigned i = 0;
i <
ndof;
i++)
5047 std::list<std::pair<unsigned, double>>*
list_pt;
5089 for (
unsigned i = 0;
i <
nvar;
i++)
5102 for (
unsigned j = 0;
j <
nvar;
j++)
5125 std::make_pair(
unknown, value));
5137 std::make_pair(eqn_number, value));
5206 value[
m] =
new double[nnz[
m]];
5228 std::list<std::pair<unsigned, double>>::iterator
it =
5301 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
5302 pause(
"Check memory usage now.");
5335 unsigned long el_lo = 0;
5378 <<
"row_or_column_start.size() "
5380 <<
"column_or_row_index.size() "
5390 <<
"value.size() " << value.
size() <<
" does not equal "
5391 <<
"column_or_row_index.size() "
5421 for (
unsigned i = 0;
i <
ndof;
i++)
5489 for (
unsigned i = 0;
i <
nvar;
i++)
5502 for (
unsigned j = 0;
j <
nvar;
j++)
5524 for (
unsigned k = 0;
k <= size;
k++)
5529 std::make_pair(
unknown, value));
5545 for (
unsigned k = 0;
k <= size;
k++)
5550 std::make_pair(eqn_number, value));
5617 for (
unsigned long i = 0;
i <
ndof;
i++)
5654 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
5655 pause(
"Check memory usage now.");
5688 unsigned long el_lo = 0;
5732 <<
"row_or_column_start.size() "
5734 <<
"column_or_row_index.size() "
5744 <<
"value.size() " << value.
size() <<
" does not equal "
5745 <<
"column_or_row_index.size() "
5778 for (
unsigned i = 0;
i <
ndof;
i++)
5847 for (
unsigned i = 0;
i <
nvar;
i++)
5860 for (
unsigned j = 0;
j <
nvar;
j++)
5881 const unsigned size =
5884 for (
unsigned k = 0;
k <= size;
k++)
5906 const unsigned size =
5908 for (
unsigned k = 0;
k <= size;
k++)
5980 for (
unsigned long i = 0;
i <
ndof;
i++)
6017 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
6018 pause(
"Check memory usage now.");
6051 unsigned long el_lo = 0;
6095 <<
"row_or_column_start.size() "
6097 <<
"column_or_row_index.size() "
6107 <<
"value.size() " << value.
size() <<
" does not equal "
6108 <<
"column_or_row_index.size() "
6141 for (
unsigned i = 0;
i <
ndof;
i++)
6225 for (
unsigned i = 0;
i <
nvar;
i++)
6238 for (
unsigned j = 0;
j <
nvar;
j++)
6254 const unsigned size =
ncoef[
m][eqn_number];
6261 [
m][eqn_number] != 0)
6287 for (
unsigned k = 0;
k <= size;
k++)
6293 [
m][eqn_number] ==
ncoef[
m][eqn_number])
6300 for (
unsigned c = 0; c <
ncoef[
m][eqn_number]; c++)
6315 unsigned entry =
ncoef[
m][eqn_number];
6335 for (
unsigned k = 0;
k <= size;
k++)
6430 for (
unsigned long i = 0;
i <
ndof;
i++)
6475 oomph_info <<
"Pausing at end of sparse assembly." << std::endl;
6476 pause(
"Check memory usage now.");
6488 const unsigned&
el_lo,
6489 const unsigned&
el_hi,
6510 for (
unsigned i = 0;
i <
nvar;
i++)
6513 unsigned global_eqn_number =
6561 <<
"This is usually a sign that the problem distribution \n"
6562 <<
"or the load balancing have gone wrong.";
6564 "Problem::parallel_sparse_assemble()",
6571 unsigned long el_lo = 0;
6607 <<
"row_or_column_start.size() " << row_start.
size()
6608 <<
" does not equal "
6619 <<
"value.size() " << values.
size() <<
" does not equal "
6756 for (
unsigned i = 0;
i <
nvar;
i++)
6759 unsigned global_eqn_number =
6766 int eqn_number =
right / 2;
6767 while (
my_eqns[eqn_number] != global_eqn_number)
6785 for (
unsigned j = 0;
j <
nvar;
j++)
6811 <<
"Internal Error: " << std::endl
6812 <<
"Could not find global equation number "
6813 << global_eqn_number
6814 <<
" in my_eqns vector of equation numbers but\n"
6815 <<
"at least one entry in the residual vector is nonzero.";
6825 if (
my_eqns[eqn_number] > global_eqn_number)
6844 for (
unsigned j = 0;
j <
nvar;
j++)
6860 const unsigned size =
ncoef[
m][eqn_number];
6867 [
m][eqn_number] != 0)
6892 for (
unsigned k = 0;
k <= size;
k++)
6898 [
m][eqn_number] ==
ncoef[
m][eqn_number])
6905 for (
unsigned c = 0; c <
ncoef[
m][eqn_number]; c++)
6920 unsigned entry =
ncoef[
m][eqn_number];
6983 oomph_info <<
"\nCPU for residual computation (loc/max/min/imbal): ";
6987 oomph_info <<
"\nCPU for Jacobian computation (loc/max/min/imbal): ";
7014 for (
unsigned c = 0; c <
ncoef[
m][
e]; c++)
7073 for (
unsigned p = 0;
p <
nproc;
p++)
7091 for (
unsigned p = 0;
p <
nproc;
p++)
7134 for (
unsigned p = 0;
p <
nproc;
p++)
7152 for (
unsigned p = 0;
p <
nproc;
p++)
7172 for (
unsigned p = 0;
p <
nproc;
p++)
7218 for (
unsigned p = 0;
p <
nproc;
p++)
7258 for (
unsigned p = 0;
p <
nproc;
p++)
7457 row_start[
m][0] = 0;
7462 for (
unsigned p = 0;
p <
nproc;
p++)
7478 row_start[
m][
i] = 0;
7482 for (
unsigned p = 0;
p <
nproc;
p++)
7537 for (
unsigned p = 0;
p <
nproc;
p++)
7562 for (
unsigned c = 0; c <
n_chunk; c++)
7629 for (
unsigned p = 0;
p <
nproc;
p++)
7644 for (
unsigned c = 0; c <
n_chunk; c++)
7651 values[
m] =
new double[nnz[
m]];
7656 for (
unsigned c = 0; c <
n_chunk; c++)
7659 for (
unsigned i = 0;
i <
nc;
i++)
7672 unsigned g = row_start[
m][0];
7673 row_start[
m][0] = 0;
7676 unsigned h = g + row_start[
m][
i];
7677 row_start[
m][
i] = g;
7691 for (
unsigned p = 0;
p <
nproc;
p++)
7706 for (
unsigned p = 0;
p <
nproc;
p++)
7717 for (
unsigned p = 0;
p <
nproc;
p++)
7733 for (
unsigned p = 0;
p <
nproc;
p++)
7784 oomph_info <<
"CPU for residual distribut. (loc/max/min/imbal): ";
7788 oomph_info <<
"CPU for Jacobian distribut. (loc/max/min/imbal): ";
7808 OomphLibWarning(
"This is unlikely to work with a distributed problem",
7809 " Problem::get_fd_jacobian()",
7824 const double FD_step = 1.0e-8;
7920 const double FD_step = 1.0e-8;
7980 for (
unsigned i = 0;
i <
n_vec;
i++)
7990 for (
unsigned i = 0;
i <
n_vec;
i++)
8025 for (
unsigned l = 0;
l <
n_var;
l++)
8028 const unsigned long eqn_number =
8032 for (
unsigned i = 0;
i <
n_vec;
i++)
8034 C_local(
i,
l) = C[
i].global_value(eqn_number);
8043 for (
unsigned l = 0;
l <
n_var;
l++)
8045 const unsigned long eqn_number =
8048 for (
unsigned i = 0;
i <
n_vec;
i++)
8084 for (
unsigned i = 0;
i <
n_vec;
i++)
8102 for (
unsigned i = 0;
i <
n_vec;
i++)
8118 for (
unsigned i = 0;
i <
n_vec;
i++)
8127 for (
unsigned i = 0;
i <
n_vec;
i++)
8141 for (
unsigned long e = 0;
e < n_element;
e++)
8160 for (
unsigned n = 0;
n <
n_var;
n++)
8167 for (
unsigned i = 0;
i <
n_vec;
i++)
8170 for (
unsigned n = 0;
n <
n_var;
n++)
8175 C_mult[
i] * C[
i].global_value(eqn_number);
8186 for (
unsigned n = 0;
n <
n_var;
n++)
8194 for (
unsigned n = 0;
n <
n_var;
n++)
8198 for (
unsigned m = 0;
m <
n_var;
m++)
8220 for (
unsigned i = 0;
i <
n_vec;
i++)
8222 product[
i].sum_all_halo_and_haloed_values();
8233 Vector<std::complex<double>>& alpha,
8300 Vector<std::complex<double>>& eigenvalue,
8365 Vector<std::complex<double>>& eigenvalue,
8431 const double&
shift)
8447 <<
"The distributions of the jacobian and mass matrix are\n"
8448 <<
"not the same and they must be.\n";
8457 <<
"mass_matrix has a distribution, but the number of rows is not "
8458 <<
"equal to the number of degrees of freedom in the problem.";
8467 <<
"main_matrix has a distribution, but the number of rows is not "
8468 <<
"equal to the number of degrees of freedom in the problem.";
8478 error_stream <<
"The distribution of the jacobian and mass matrix must "
8479 <<
"both be setup or both not setup";
8636 (*Saved_dof_pt)[
i] = *(this->
Dof_pt[
i]);
8650 for (
unsigned long n = 0;
n <
n_dof;
n++)
8652 (*Saved_dof_pt)[
n] =
dof(
n);
8666 "There are no stored values, use store_current_dof_values()\n",
8681 throw OomphLibError(
"The number of stored values is not equal to the "
8682 "current number of dofs\n",
8702 throw OomphLibError(
"The number of stored values is not equal to the "
8703 "current number of dofs\n",
8709 for (
unsigned long n = 0;
n <
n_dof;
n++)
8711 dof(
n) = (*Saved_dof_pt)[
n];
8729 std::ostringstream error_message;
8730 error_message <<
"Eigenvector has size " <<
eigenvector.nrow()
8731 <<
", not equal to the number of dofs in the problem,"
8732 <<
n_dof << std::endl;
8767 std::ostringstream error_message;
8768 error_message <<
"Eigenvector has size " <<
eigenvector.nrow()
8769 <<
", not equal to the number of dofs in the problem,"
8770 <<
n_dof << std::endl;
8839 error_stream <<
"Globally convergent Newton method has not been "
8840 <<
"implemented in parallel yet!" << std::endl;
8871 oomph_info << std::endl << std::endl << std::endl;
8872 oomph_info <<
"This is a bit bizarre: The problem has no dofs."
8875 <<
"I'll just return from the Newton solver without doing anything."
8884 oomph_info <<
"I hope this is what you intended me to do..."
8887 <<
"Note: All actions_...() functions were called"
8889 oomph_info << std::endl <<
" before returning." << std::endl;
8890 oomph_info << std::endl << std::endl << std::endl;
8926 double maxres =
dx.max();
8940 oomph_info <<
"\nInitial Maximum residuals " << maxres << std::endl;
8955 <<
"Linear problem -- convergence in one iteration assumed."
8973 oomph_info <<
"Not recomputing Jacobian! " << std::endl;
8993 oomph_info <<
"Enabling resolve" << std::endl;
9017 double*
dx_pt =
dx.values_pt();
9064 double maxres = 0.0;
9085 << maxres << std::endl
9106 <<
") has been exceeded in Newton solver." << std::endl;
9110 oomph_info <<
"Reached max. number of iterations ("
9146 oomph_info <<
"Time outside linear solver : "
9149 <<
" %" << std::endl;
9156 oomph_info <<
"Time outside linear solver : "
9157 <<
"[too fast]" << std::endl;
9181 for (
unsigned i = 0;
i <
n_dof;
i++)
9188 for (
unsigned i = 0;
i <
n_dof;
i++)
9194 for (
unsigned i = 0;
i <
n_dof;
i++)
9201 warn_message <<
"WARNING: Non-negative slope, probably due to a "
9202 <<
" roundoff \nproblem in the linesearch: slope=" <<
slope
9205 "Problem::globally_convergent_line_search()",
9209 for (
unsigned i = 0;
i <
n_dof;
i++)
9216 double lambda = 1.0;
9219 for (
unsigned i = 0;
i <
n_dof;
i++)
9228 for (
unsigned i = 0;
i <
n_dof;
i++)
9239 warn_message <<
"WARNING: Line search converged on x only!\n";
9241 "Problem::globally_convergent_line_search()",
9248 oomph_info <<
"Returning from linesearch with lambda=" << lambda
9344 <<
"USER-DEFINED ERROR IN NEWTON SOLVER " << std::endl;
9346 if (
error.linear_solver_error())
9348 oomph_info <<
"ERROR IN THE LINEAR SOLVER" << std::endl;
9354 <<
") REACHED WITHOUT CONVERGENCE " << std::endl;
9366 error_stream <<
"Error occured in Newton solver. " << std::endl;
9471 double maxres = y.
max();
9512 oomph_info <<
"Initial Maximum residuals " << maxres << std::endl;
9542 ->solve_for_two_rhs(
this, y,
input_z, z);
9649 double maxres = y.
max();
9889 if (n_sub_mesh == 0)
9893 if (spine_mesh_pt->does_pointer_correspond_to_spine_data(
parameter_pt))
9908 if (spine_mesh_pt->does_pointer_correspond_to_spine_data(
10043#ifdef OOMPH_HAS_MPI
10083 <<
"Warning: This function is called after spatially adapting the\n"
10084 <<
"eigenfunction associated with a pitchfork bifurcation and should\n"
10085 <<
"ensure that the exact (anti-)symmetries of problem are enforced\n"
10086 <<
"within that eigenfunction. It is problem specific and must be\n"
10087 <<
"filled in by the user if required.\n"
10088 <<
"A sign of problems is if the slack paramter gets too large and\n"
10089 <<
"if the solution at the Pitchfork is not symmetric.\n";
10092 "Problem::symmetrise_eigenfunction_for_adaptive_pitchfork_tracking()",
10267 const double& omega,
10322 std::ostringstream error_message;
10324 <<
"The parameter addressed by " <<
parameter_pt <<
" with the value "
10326 <<
"\n is supposed to be used for arc-length contiunation,\n"
10327 <<
" but it is stored in a Data object used by the problem.\n\n"
10328 <<
"This is bad for two reasons:\n"
10329 <<
"1. If it's a variable in the problem, it must already have an\n"
10330 "associated equation, so it can't be used for continuation;\n"
10331 <<
"2. The problem data will be reorganised in memory during "
10333 <<
" which means that the pointer will become invalid.\n\n"
10334 <<
"If you are sure that this is what you want to do you must:\n"
10335 <<
"A. Ensure that the value is pinned (don't worry we'll shout again "
10337 <<
"B. Use the alternative interface\n"
10338 <<
" Problem::arc_length_step_solve(Data*,unsigned,...)\n"
10339 <<
" which uses a pointer to the data object and not the raw double "
10353 for (
unsigned i = 0;
i < n_time_steppers;
i++)
10357 continuation_time_stepper_added =
true;
10365 oomph_info <<
"Adding the continuation time stepper\n";
10381 <<
" equation numbers allocated for continuation\n";
10412 <<
" in the data object to be used for continuation\n"
10413 <<
"is not pinned, which means that it is already a\n"
10414 <<
"variable in the problem "
10415 <<
"and cannot be used for continuation.\n\n"
10416 <<
"Please correct your formulation by either:\n"
10417 <<
"A. Pinning the value"
10419 <<
"B. Using a different parameter for continuation"
10432 for (
unsigned i = 0;
i < n_time_steppers;
i++)
10436 continuation_time_stepper_added =
true;
10444 oomph_info <<
"Adding the continuation time stepper\n";
10461 <<
" equation numbers allocated for continuation\n";
10495 if (n_sub_mesh == 0)
10499 spine_mesh_pt->set_consistent_pinned_spine_values_for_continuation(
10513 spine_mesh_pt->set_consistent_pinned_spine_values_for_continuation(
10590#ifdef OOMPH_HAS_MPI
10616#ifdef OOMPH_HAS_MPI
10617 <<
", in total (over all processors).\n";
10626 oomph_info <<
"\n \n Solution is fully converged in "
10627 <<
"Problem::newton_solver(). \n \n ";
10673 unsigned count = 0;
10691 std::ostringstream error_message;
10692 error_message <<
"DESIRED ARC-LENGTH STEP " <<
Ds_current
10693 <<
" HAS FALLEN BELOW MINIMUM TOLERANCE, " <<
Minimum_ds
10727 if (
error.linear_solver_error())
10731 <<
"USER-DEFINED ERROR IN NEWTON SOLVER " << std::endl;
10732 oomph_info <<
"ERROR IN THE LINEAR SOLVER" << std::endl;
10740 oomph_info <<
"STEP REJECTED DUE TO NEWTON SOLVER --- TRYING AGAIN"
10750 <<
"STEP REJECTED DUE TO INVERTED ELEMENTS --- TRYING AGAIN"
10775 std::string error_message =
10776 "The sign of the jacobian is zero after a linear solve\n";
10777 error_message +=
"Either the matrix is singular (unlikely),\n";
10778 error_message +=
"or the linear solver cannot compute the "
10779 "determinant of the matrix;\n";
10780 error_message +=
"e.g. an iterative linear solver.\n";
10782 "If the latter, bifurcation detection must be via an eigensolver\n";
10784 "Problem::arc_length_step_solve",
10819 <<
"-----------------------------------------------------------";
10821 <<
"SIGN CHANGE IN DETERMINANT OF JACOBIAN: " << std::endl;
10822 message <<
"BIFURCATION OR TURNING POINT DETECTED BETWEEN "
10828 <<
"-----------------------------------------------------------"
10836 std::ios_base::app);
10953 throw OomphLibError(
"Explicit time stepper pointer is null in problem.",
11039 <<
"USER-DEFINED ERROR IN NEWTON SOLVER " << std::endl;
11041 if (
error.linear_solver_error())
11043 oomph_info <<
"ERROR IN THE LINEAR SOLVER" << std::endl;
11049 <<
") REACHED WITHOUT CONVERGENCE " << std::endl;
11060 error_stream <<
"Error occured in unsteady Newton solver. " << std::endl;
11087 const double& epsilon)
11107 const double& epsilon,
11129 unsigned adaptive_flag = 0;
11196 if (
error.linear_solver_error() ||
11199 std::string error_message =
"USER-DEFINED ERROR IN NEWTON SOLVER\n";
11200 error_message +=
"ERROR IN THE LINEAR SOLVER\n";
11209 oomph_info <<
"TIMESTEP REJECTED DUE TO THE NEWTON SOLVER"
11220 oomph_info <<
"TIMESTEP REJECTED DUE TO INVERTED ELEMENTS" << std::endl;
11265 <<
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n"
11266 <<
"Estimated timestepping error is " <<
error <<
"\n"
11271 if (
error > epsilon)
11274 <<
" exceeds tolerance " << epsilon <<
"\n";
11282 oomph_info <<
" ...but we're not rejecting the timestep\n";
11285 <<
"Note: This behaviour can be adjusted by changing the\n"
11286 <<
"protected boolean\n"
11287 <<
" Problem::Keep_temporal_error_below_tolerance\n\n"
11288 <<
"Also, if you are noticing that many of your timesteps result\n"
11289 <<
"in error > tolerance, try reducing the target error with\n"
11290 <<
"respect to the error tolerance by reducing the value of\n"
11291 <<
"Target_error_safety_factor from its default value of 1.0\n"
11292 <<
"using the access function\n"
11293 <<
" target_error_safety_factor() = 0.5 (e.g.)\n"
11294 <<
"The default strategy (Target_error_safety_factor=1.0) tries\n"
11295 <<
"to suggest a timestep which will produce an error equal to\n"
11296 <<
"the error tolerance `epsilon` which risks error > tolerance\n"
11297 <<
"quite often. Setting the safety factor to too small a value\n"
11298 <<
"will make the timesteps unnecessarily small; too large will\n"
11299 <<
"not address the issue -- neither is optimal and a problem\n"
11300 <<
"dependent compromise is needed.\n"
11301 <<
"for more info see:\n"
11302 <<
" Mayr et al. (2018), p5,9, DOI:10.1016/j.finel.2017.12.002\n"
11303 <<
" Harrier et al. (1993), p168, ISBN:978-3-540-56670-0\n"
11304 <<
" Söderlind (2002), (2.7) on p5, DOI:10.1023/A:1021160023092\n";
11307 <<
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n"
11324 oomph_info <<
"Tried to increase dt by the ratio "
11327 <<
"). Attempting to increase by the maximum ratio instead."
11340 <<
"Warning: Adaptation of timestep to ensure satisfaction\n"
11341 <<
" of error bounds during adaptive timestepping\n"
11342 <<
" would lower dt below \n"
11343 <<
" Problem::Minimum_dt_but_still_proceed="
11345 <<
" ---> We're continuing with present timestep.\n"
11356 <<
" which is less than the minimum scaling factor "
11358 oomph_info <<
"TIMESTEP REJECTED" << std::endl;
11367 <<
" which is above the maximum (" <<
Maximum_dt
11368 <<
"). I increased it to the maximum value instead.";
11373 std::ostringstream
err;
11375 <<
" which is less than the minimum dt (" <<
Minimum_dt <<
")."
11397 for (
unsigned i = 0;
i <
ni;
i++)
11402#ifdef OOMPH_HAS_MPI
11417 while (reject_timestep);
11444 const double& epsilon,
11464 oomph_info <<
"No spatial refinement allowed; max_adapt=0\n";
11479#ifdef OOMPH_HAS_MPI
11501 oomph_info <<
"Mesh was adapted but re-solve has been suppressed."
11507 <<
"Mesh was adapted --> we'll re-solve for current timestep."
11516 bool shift =
false;
11525 oomph_info <<
"Re-assigning initial condition at time="
11548 oomph_info <<
"Mesh wasn't adapted --> we'll accept spatial refinement."
11648 if (n_sub_mesh == 0)
11675#ifdef OOMPH_HAS_MPI
11679 warning_stream <<
"This has not been comprehensively tested for "
11680 "distributed problems.\n"
11681 <<
"I'm sure that I need to worry about external halo and "
11682 "external elements."
11732 std::string
err =
"Prediction by explicit step only works for "
11733 "problems with a simple time";
11734 err +=
"stepper. I think implementing anything more general will";
11735 err +=
"require a rewrite of explicit time steppers. - David";
11756 std::string
err =
"Requested predictions by explicit step but explicit";
11757 err +=
" predictor pt is null.";
11765 throw OomphLibError(
"Problem has explicit time stepper other than "
11766 "predictor, not sure how to handle this yet ??ds",
11798 using namespace StringConversion;
11799 std::string
err =
"Predictor landed at the wrong time!";
11834#ifdef OOMPH_HAS_MPI
11837 throw OomphLibError(
"Not yet implemented for distributed problems",
11848 std::string
err =
"Not implemented for multiple time steppers";
11858 for (
unsigned i = 0;
i <
ndof();
i++)
11886 elem_pt->enable_mass_matrix_reuse();
11911 elem_pt->disable_mass_matrix_reuse();
11939 oomph_info <<
"Copying an unsteady problem." << std::endl;
11945 for (
unsigned i = 0;
i <
n_dt;
i++)
11965 if (nmesh == 0) nmesh = 1;
11966 for (
unsigned m = 0;
m < nmesh;
m++)
11975 std::ostringstream error_message;
11976 error_message <<
"Number of nodes in copy " <<
n_node
11977 <<
" not equal to the number in the original "
11986 for (
unsigned long i = 0;
i <
n_node;
i++)
12015 std::ostringstream error_message;
12016 error_message <<
"Number of global data in copy " <<
n_global
12017 <<
" not equal to the number in the original "
12034 for (
unsigned m = 0;
m < nmesh;
m++)
12049 std::ostringstream error_message;
12050 error_message <<
"Number of internal data in copy " <<
n_internal
12051 <<
" not equal to the number in the original "
12078 <<
"This function must be overloaded in your specific problem, and must\n"
12079 <<
"create an exact copy of your problem. Usually this will be achieved\n"
12080 <<
"by a call to the constructor with exactly the same arguments as "
12097 dump_file << std::max(
unsigned(1),
n_mesh) <<
" # number of (sub)meshes "
12109 <<
" # uniform refinement when pruned " << std::endl;
12113 dump_file << 0 <<
" # (fake) uniform refinement when pruned "
12116 dump_file << 9999 <<
" # test flag for end of sub-meshes " << std::endl;
12130 <<
" # uniform refinement when pruned " << std::endl;
12134 dump_file << 0 <<
" # (fake) uniform refinement when pruned "
12138 dump_file << 9999 <<
" # test flag for end of sub-meshes " << std::endl;
12141#ifdef OOMPH_HAS_MPI
12149 for (
unsigned e = 0;
e <
n;
e++)
12186 dump_file <<
n <<
" # Number of base elements; partitioning follows.\n";
12187 for (
unsigned e = 0;
e <
n;
e++)
12191 dump_file <<
"8888 #test flag for end of base element distribution\n";
12205#ifdef OOMPH_HAS_TRIANGLE_LIB
12210#ifdef OOMPH_HAS_MPI
12213 if (
mmesh_pt->is_mesh_distributed())
12237#ifdef OOMPH_HAS_TRIANGLE_LIB
12243#ifdef OOMPH_HAS_MPI
12246 if (
mmesh_pt->is_mesh_distributed())
12273 dump_file <<
n_dt <<
" # Number of timesteps " << std::endl;
12274 for (
unsigned i = 0;
i <
n_dt;
i++)
12283 dump_file <<
"0.0 # Dummy time from steady run " << std::endl;
12285 dump_file <<
"0 # Dummy number of timesteps from steady run" << std::endl;
12290 if (nmesh == 0) nmesh = 1;
12291 for (
unsigned m = 0;
m < nmesh;
m++)
12324 warn_message <<
"Restart file isn't open -- I'm assuming that this is\n";
12325 warn_message <<
"because we're restarting on a larger number of\n";
12326 warn_message <<
"processor than were in use when the restart data was \n";
12353 std::ostringstream error_message;
12355 <<
"Number of sub-meshes specified in restart file, "
12356 <<
n_submesh_read <<
" doesn't \n match the my number of sub-meshes,"
12358 <<
"Make sure all sub-meshes have been added to the global mesh\n"
12359 <<
"when calling the Problem::dump() function.\n";
12373#ifdef OOMPH_HAS_MPI
12398 error_stream <<
"Nonzero uniform-refinement-when-pruned specified\n"
12399 <<
"even though mesh is not tree-based. Odd. May want\n"
12400 <<
"to check this carefully before disabling this \n"
12401 <<
"warning/error -- most likely if/when we start to\n"
12402 <<
"prune unstructured meshes [though I can't see why\n"
12403 <<
"we would want to do this, given that they are \n"
12404 <<
"currently totally re-generated...]\n";
12420#ifdef OOMPH_HAS_MPI
12452#ifdef OOMPH_HAS_MPI
12463#ifdef OOMPH_HAS_MPI
12504 std::ostringstream error_message;
12505 error_message <<
"Number of uniform refinements was not consistent \n"
12506 <<
"for following meshes during restart on processor \n"
12507 <<
"on which restart file could be opened:\n";
12514 error_message <<
"Sub-mesh: " <<
i <<
"; local nrefinement: "
12516 <<
"; global/synced nrefinement: "
12545 std::ostringstream error_message;
12547 <<
"Error in reading restart data: Uniform refinement when pruned \n"
12548 <<
"flags should be followed by 9999.\n";
12562#ifdef OOMPH_HAS_MPI
12618 <<
"The number of base elements in the mesh is 0,\n"
12619 <<
" but the restart file indicates that there are "
12621 <<
"This could be because the restart file was \n"
12622 <<
"generated by using code without MPI.\n"
12624 <<
"The best fix is to include two additional lines\n"
12625 <<
"in the restart file: \n\n"
12626 <<
"0 # Number of base elements; partitioning follows.\n"
12627 <<
"8888 # Test flag for end of base element distribution\n"
12629 <<
"These lines go after the flag 9999 that indicates\n"
12630 <<
"the end of the submesh information.\n"
12632 <<
"The file will now continue to be read assuming that\n"
12633 <<
"the base element information is not present.\n"
12634 <<
"If you get strange results then please look carefully\n"
12635 <<
"at the restart file. The safest thing to do is to \n"
12636 <<
"ensure that the restart file was generated by code\n"
12637 <<
"compiled and run with the same parallel options.\n";
12649 std::ostringstream error_message;
12651 <<
" base elements \n"
12652 <<
"though we only have " <<
nbase
12653 <<
" base elements in mesh.\n";
12666 for (
unsigned e = 0;
e <
nbase;
e++)
12690 std::ostringstream error_message;
12692 <<
"Error in reading restart data: Target proc for base elements \n"
12693 <<
"should be followed by 8888.\n";
12714 for (
unsigned e = 0;
e <
nel;
e++)
12743 "el_number_in_base_mesh_plus_one=0 for bulk",
12789 oomph_info <<
"Doing load balancing after pruning\n";
12795 oomph_info <<
"Done load balancing after pruning\n";
12799 oomph_info <<
"No need for load balancing after pruning\n";
12847#ifdef OOMPH_HAS_TRIANGLE_LIB
12852#ifdef OOMPH_HAS_MPI
12855 if (
mmesh_pt->is_mesh_distributed())
12866#ifdef OOMPH_HAS_MPI
12870 if (
mmesh_pt->is_mesh_distributed())
12872 mmesh_pt->reestablish_distribution_info_for_restart(
12903#ifdef OOMPH_HAS_TRIANGLE_LIB
12909#ifdef OOMPH_HAS_MPI
12912 if (
mmesh_pt->is_mesh_distributed())
12924#ifdef OOMPH_HAS_MPI
12928 if (
mmesh_pt->is_mesh_distributed())
12930 mmesh_pt->reestablish_distribution_info_for_restart(
12965 <<
"I've just read in some unstructured meshes and have, in\n"
12966 <<
"the process, totally re-generated their nodes and elements.\n"
12967 <<
"This may create dangling pointers that still point to the\n"
12968 <<
"old nodes and elements, e.g. because FaceElements were\n"
12969 <<
"attached to these meshes or pointers to nodes and elements\n"
12970 <<
"were stored somewhere. FaceElements should therefore be\n"
12971 <<
"removed before reading in these meshes, using an overloaded\n"
12972 <<
"version of the function\n\n"
12973 <<
" Problem::actions_before_read_unstructured_meshes()\n\n"
12974 <<
"and then re-attached using an overloaded version of\n\n"
12975 <<
" Problem::actions_after_read_unstructured_meshes().\n\n"
12976 <<
"The required content of these functions is likely to be "
12978 <<
"to the Problem::actions_before_adapt() and \n"
12979 <<
"Problem::actions_after_adapt() that would be required in\n"
12980 <<
"a spatially adaptive computation. If these functions already\n"
12981 <<
"exist and perform the required actions, the \n"
12982 <<
"actions_before/after_read_unstructured_meshes() functions\n"
12983 <<
"can remain empty because the former are called automatically.\n"
12984 <<
"In this case, this warning my be suppressed by setting the\n"
12985 <<
"public boolean\n\n"
12987 "Problem::Suppress_warning_about_actions_before_read_"
12988 "unstructured_meshes\n\n"
12989 <<
"to true." << std::endl;
12998 oomph_info <<
"\nNumber of equations in Problem::read(): "
13006#ifdef OOMPH_HAS_MPI
13013 oomph_info <<
"Restart file exists" << std::endl;
13014#ifdef OOMPH_HAS_MPI
13061 oomph_info <<
"Restart file does not exist" << std::endl;
13062#ifdef OOMPH_HAS_MPI
13074#ifdef OOMPH_HAS_MPI
13090#ifdef OOMPH_HAS_MPI
13096 std::ostringstream error_message;
13097 error_message <<
"Synchronisation of temporal restart data \n"
13098 <<
"required even though Problem hasn't been distributed "
13159 std::ostringstream error_message;
13161 <<
"Synchronisation of temporal restart data \n"
13162 <<
"required even though we don't have mpi support -- very odd!\n";
13191 if (nmesh == 0) nmesh = 1;
13192 for (
unsigned m = 0;
m < nmesh;
m++)
13238#ifdef OOMPH_HAS_TRIANGLE_LIB
13248 mmesh_pt->update_polyline_representation_from_restart();
13273 std::ostringstream error_message;
13274 error_message <<
"The number of global data " <<
Nglobal
13275 <<
" is not equal to that specified in the input file "
13354 <<
"\n ERROR: Failed Mesh::self_test() for single mesh in problem"
13366 oomph_info <<
"\n ERROR: Failed Mesh::self_test() for mesh imesh"
13367 <<
imesh << std::endl;
13381 <<
"\n ERROR: Failed Data::self_test() for global data iglobal"
13387#ifdef OOMPH_HAS_MPI
13419 const unsigned& bifurcation_type,
13431 double omega = 0.0;
13433 if (bifurcation_type == 3)
13439 double sigma = 0.0;
13440 if (bifurcation_type == 2)
13442 sigma = this->
dof(this->
ndof() - 1);
13456 for (
unsigned c = 0; c <
n_copies; c++)
13477 if (
mmesh_pt->is_adaptation_enabled())
13484 mmesh_pt->refine_base_mesh_as_in_reference_mesh(
13490 <<
"Info/Warning: Mesh in orginal problem is not refineable."
13496 oomph_info <<
"Info/Warning: Mesh adaptation is disabled in copy."
13502 oomph_info <<
"Info/Warning: Mesh cannot be adapted in copy."
13517 if (
mmesh_pt->is_adaptation_enabled())
13524 mmesh_pt->refine_base_mesh_as_in_reference_mesh(
13529 oomph_info <<
"Info/Warning: Mesh in orginal problem is not "
13537 <<
"Info/Warning: Mesh adaptation is disabled in copy."
13543 oomph_info <<
"Info/Warning: Mesh cannot be adapted in copy."
13562 for (
unsigned c = 0; c <
n_copies; c++)
13570 error_stream <<
"Number of unknowns in the problem copy " << c <<
" "
13571 <<
"not equal to number in the original:\n"
13572 << this->
ndof() <<
" (original) "
13587 if (bifurcation_type == 2)
13590 ->symmetrise_eigenfunction_for_adaptive_pitchfork_tracking();
13598 for (
unsigned c = 0; c <
n_copies; c++)
13610 error_stream <<
"Problems do not have the same number of meshes\n"
13628 <<
" does not have the same number of elements in the "
13652 for (
unsigned c = 0; c <
n_copies; c++)
13658 if (bifurcation_type == 2)
13661 ->symmetrise_eigenfunction_for_adaptive_pitchfork_tracking();
13665 for (
unsigned c = 0; c <
n_copies; c++)
13672 switch (bifurcation_type)
13683 this->
dof(this->
ndof() - 1) = sigma;
13694 error_stream <<
"Bifurcation type " << bifurcation_type
13696 <<
"1: Fold, 2: Pitchfork, 3: Hopf\n";
13753 if (bifurcation_type != 0)
13765 for (
unsigned c = 0; c < 2; c++)
13788 if (
mmesh_pt->is_adaptation_enabled())
13798 <<
"Info/Warning: Adaptive Continuation is broken in "
13799 <<
"SolidElement" << std::endl;
13801 mmesh_pt->refine_base_mesh_as_in_reference_mesh(
13806 oomph_info <<
"Info/Warning: Mesh in orginal problem is not "
13814 <<
"Info/Warning: Mesh adaptation is disabled in copy."
13828 <<
"Info/Warning: Adaptive Continuation is broken in "
13829 <<
"SolidElement" << std::endl;
13835 std::ofstream
tri_dump(
"triangle_mesh.dmp");
13838 std::ifstream
tri_read(
"triangle_mesh.dmp");
13852 for (
unsigned i = 0;
i <
n_dim; ++
i)
13861 <<
"Info/warning: Original Mesh is not TriangleBased\n"
13862 <<
"... but the copy is!" << std::endl;
13867 oomph_info <<
"Info/Warning: Mesh cannot be adapted in copy."
13882 if (
mmesh_pt->is_adaptation_enabled())
13892 <<
"Info/Warning: Adaptive Continuation is broken in "
13893 <<
"SolidElement" << std::endl;
13896 mmesh_pt->refine_base_mesh_as_in_reference_mesh(
13901 oomph_info <<
"Info/Warning: Mesh in orginal problem is "
13909 <<
"Info/Warning: Mesh adaptation is disabled in copy."
13923 <<
"Info/Warning: Adaptive Continuation is broken in "
13924 <<
"SolidElement" << std::endl;
13930 std::ofstream
tri_dump(
"triangle_mesh.dmp");
13933 std::ifstream
tri_read(
"triangle_mesh.dmp");
13946 for (
unsigned i = 0;
i <
n_dim; ++
i)
13955 <<
"Info/warning: Original Mesh is not TriangleBased\n"
13956 <<
"... but the copy is!" << std::endl;
13961 oomph_info <<
"Info/Warning: Mesh cannot be adapted in copy."
13988 error_stream <<
"Number of unknowns in the problem copy " << c <<
" "
13989 <<
"not equal to number in the original:\n"
13990 << this->
ndof() <<
" (original) "
14041 oomph_info <<
"Adapting problem:" << std::endl;
14042 oomph_info <<
"=================" << std::endl;
14053 double t_end = 0.0;
14077 if (
mmesh_pt->is_adaptation_enabled())
14083 mmesh_pt->spatial_error_estimator_pt();
14108 mmesh_pt->max_error() = std::fabs(*std::max_element(
14111 mmesh_pt->min_error() = std::fabs(*std::min_element(
14115 <<
mmesh_pt->min_error() << std::endl
14137 oomph_info <<
"Time for complete mesh adaptation "
14138 <<
"(but excluding comp of error estimate): "
14145 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14151 oomph_info <<
"Info/Warning: Mesh cannot be adapted" << std::endl;
14169 mmesh_pt->spatial_error_estimator_pt();
14180 if (
mmesh_pt->is_adaptation_enabled())
14210 <<
mmesh_pt->min_error() << std::endl;
14232 oomph_info <<
"Time for complete mesh adaptation "
14233 <<
"(but excluding comp of error estimate): "
14240 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14246 oomph_info <<
"Info/Warning: Mesh cannot be adapted." << std::endl;
14259 oomph_info <<
"Total time for actual adaptation "
14260 <<
"(all meshes; incl error estimates): " <<
t_end -
t_start
14276 oomph_info <<
"About to start re-assigning eqn numbers "
14277 <<
"with Problem::assign_eqn_numbers() at end of "
14278 <<
"Problem::adapt().\n";
14289 oomph_info <<
"Time for re-assigning eqn numbers with "
14290 <<
"Problem::assign_eqn_numbers() at end of Problem::adapt(): "
14318 if (bifurcation_type != 0)
14326 oomph_info <<
"p-adapting problem:" << std::endl;
14327 oomph_info <<
"===================" << std::endl;
14338 double t_end = 0.0;
14362 if (
mmesh_pt->is_p_adaptation_enabled())
14368 mmesh_pt->spatial_error_estimator_pt();
14393 mmesh_pt->max_error() = std::fabs(*std::max_element(
14396 mmesh_pt->min_error() = std::fabs(*std::min_element(
14400 <<
mmesh_pt->min_error() << std::endl
14422 oomph_info <<
"Time for complete mesh adaptation "
14423 <<
"(but excluding comp of error estimate): "
14430 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14436 oomph_info <<
"Info/Warning: Mesh cannot be adapted" << std::endl;
14454 mmesh_pt->spatial_error_estimator_pt();
14465 if (
mmesh_pt->is_p_adaptation_enabled())
14495 <<
mmesh_pt->min_error() << std::endl;
14517 oomph_info <<
"Time for complete mesh adaptation "
14518 <<
"(but excluding comp of error estimate): "
14525 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14531 oomph_info <<
"Info/Warning: Mesh cannot be adapted." << std::endl;
14544 oomph_info <<
"Total time for actual adaptation "
14545 <<
"(all meshes; incl error estimates): " <<
t_end -
t_start
14561 oomph_info <<
"About to start re-assigning eqn numbers "
14562 <<
"with Problem::assign_eqn_numbers() at end of "
14563 <<
"Problem::adapt().\n";
14574 oomph_info <<
"Time for re-assigning eqn numbers with "
14575 <<
"Problem::assign_eqn_numbers() at end of Problem::adapt(): "
14597 oomph_info <<
"Adapting problem:" << std::endl;
14598 oomph_info <<
"=================" << std::endl;
14618 if (
mmesh_pt->is_adaptation_enabled())
14629 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14635 oomph_info <<
"Info/Warning: Mesh cannot be adapted" << std::endl;
14650 if (
mmesh_pt->is_adaptation_enabled())
14661 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14667 oomph_info <<
"Info/Warning: Mesh cannot be adapted." << std::endl;
14706 if (
mmesh_pt->is_adaptation_enabled())
14710 mmesh_pt->spatial_error_estimator_pt();
14748 <<
mmesh_pt->min_error() << std::endl;
14752 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14758 oomph_info <<
"Info/Warning: Mesh cannot be adapted" << std::endl;
14778 mmesh_pt->spatial_error_estimator_pt();
14789 if (
mmesh_pt->is_adaptation_enabled())
14817 <<
mmesh_pt->min_error() << std::endl;
14821 oomph_info <<
"Info/Warning: Mesh adaptation is disabled."
14827 oomph_info <<
"Info/Warning: Mesh cannot be adapted." << std::endl;
14842 if (bifurcation_type != 0)
14862 mmesh_pt->spatial_error_estimator_pt();
14886 mmesh_pt->max_error() = std::fabs(*std::max_element(
14889 mmesh_pt->min_error() = std::fabs(*std::min_element(
14893 <<
mmesh_pt->min_error() << std::endl;
14910 mmesh_pt->spatial_error_estimator_pt();
14949 <<
mmesh_pt->min_error() << std::endl;
14980 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
14986 std::ostringstream error_message;
14987 error_message <<
"Problem::refine_selected_elements(...) only works for\n"
14988 <<
"multiple-mesh problems if you specify the mesh\n"
14989 <<
"number in the function argument before the Vector,\n"
14990 <<
"or a Vector of Vectors for each submesh.\n"
15026 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15032 std::ostringstream error_message;
15033 error_message <<
"Problem::refine_selected_elements(...) only works for\n"
15034 <<
"multiple-mesh problems if you specify the mesh\n"
15035 <<
"number in the function argument before the Vector,\n"
15036 <<
"or a Vector of Vectors for each submesh.\n"
15063 std::ostringstream error_message;
15064 error_message <<
"Problem only has " <<
n_mesh
15065 <<
" submeshes. Cannot refine submesh " <<
i_mesh
15079 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15111 std::ostringstream error_message;
15112 error_message <<
"Problem only has " <<
n_mesh
15113 <<
" submeshes. Cannot refine submesh " <<
i_mesh
15127 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15166 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15203 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15241 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15247 std::ostringstream error_message;
15249 <<
"Problem::p_refine_selected_elements(...) only works for\n"
15250 <<
"multiple-mesh problems if you specify the mesh\n"
15251 <<
"number in the function argument before the Vector,\n"
15252 <<
"or a Vector of Vectors for each submesh.\n"
15288 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15294 std::ostringstream error_message;
15296 <<
"Problem::p_refine_selected_elements(...) only works for\n"
15297 <<
"multiple-mesh problems if you specify the mesh\n"
15298 <<
"number in the function argument before the Vector,\n"
15299 <<
"or a Vector of Vectors for each submesh.\n"
15320 "p-refinement for multiple submeshes has not yet been tested.",
15321 "Problem::p_refine_selected_elements()",
15331 std::ostringstream error_message;
15332 error_message <<
"Problem only has " <<
n_mesh
15333 <<
" submeshes. Cannot p-refine submesh " <<
i_mesh
15347 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15373 "p-refinement for multiple submeshes has not yet been tested.",
15374 "Problem::p_refine_selected_elements()",
15384 std::ostringstream error_message;
15385 error_message <<
"Problem only has " <<
n_mesh
15386 <<
" submeshes. Cannot p-refine submesh " <<
i_mesh
15400 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15425 "p-refinement for multiple submeshes has not yet been tested.",
15426 "Problem::p_refine_selected_elements()",
15444 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15467 "p-refinement for multiple submeshes has not yet been tested.",
15468 "Problem::p_refine_selected_elements()",
15486 oomph_info <<
"Info/Warning: Mesh cannot be refined " << std::endl;
15520 double t_end = 0.0;
15525 <<
"Time for actions before adapt in Problem::refine_uniformly_aux(): "
15541 for (
unsigned i = 0;
i <
nref;
i++)
15543 mmesh_pt->refine_uniformly(doc_info);
15548 oomph_info <<
"Info/Warning: Mesh cannot be refined uniformly "
15563 for (
unsigned i = 0;
i <
nref;
i++)
15565 mmesh_pt->refine_uniformly(doc_info);
15581 oomph_info <<
"Time for mesh-level mesh refinement in "
15582 <<
"Problem::refine_uniformly_aux(): " <<
t_end -
t_start
15595 <<
"Time for actions after adapt Problem::refine_uniformly_aux(): "
15601#ifdef OOMPH_HAS_MPI
15614 oomph_info <<
"Time for Problem::prune_halo_elements_and_nodes() in "
15615 <<
"Problem::refine_uniformly_aux(): " <<
t_end -
t_start
15623 std::ostringstream error_message;
15625 <<
"Requested pruning in serial build. Ignoring the request.\n";
15627 "Problem::refine_uniformly_aux()",
15634 <<
"Number of equations after Problem::refine_uniformly_aux(): "
15640 oomph_info <<
"Time for Problem::assign_eqn_numbers() in "
15641 <<
"Problem::refine_uniformly_aux(): " <<
t_end -
t_start
15667 double t_end = 0.0;
15671 oomph_info <<
"Time for actions before adapt in "
15672 "Problem::p_refine_uniformly_aux(): "
15688 for (
unsigned i = 0;
i <
nref;
i++)
15690 mmesh_pt->p_refine_uniformly(doc_info);
15695 oomph_info <<
"Info/Warning: Mesh cannot be p-refined uniformly "
15703 "p-refinement for multiple submeshes has not yet been tested.",
15704 "Problem::p_refine_uniformly_aux()",
15715 for (
unsigned i = 0;
i <
nref;
i++)
15717 mmesh_pt->p_refine_uniformly(doc_info);
15733 oomph_info <<
"Time for mesh-level mesh refinement in "
15734 <<
"Problem::p_refine_uniformly_aux(): " <<
t_end -
t_start
15747 <<
"Time for actions after adapt Problem::p_refine_uniformly_aux(): "
15753#ifdef OOMPH_HAS_MPI
15766 oomph_info <<
"Time for Problem::prune_halo_elements_and_nodes() in "
15767 <<
"Problem::p_refine_uniformly_aux(): " <<
t_end -
t_start
15775 std::ostringstream error_message;
15777 <<
"Requested pruning in serial build. Ignoring the request.\n";
15779 "Problem::p_refine_uniformly_aux()",
15786 <<
"Number of equations after Problem::p_refine_uniformly_aux(): "
15792 oomph_info <<
"Time for Problem::assign_eqn_numbers() in "
15793 <<
"Problem::p_refine_uniformly_aux(): " <<
t_end -
t_start
15811 std::ostringstream error_message;
15812 error_message <<
"imesh " <<
i_mesh
15813 <<
" is greater than the number of sub meshes "
15825 mmesh_pt->refine_uniformly(doc_info);
15829 oomph_info <<
"Info/Warning: Mesh cannot be refined uniformly "
15855 std::ostringstream error_message;
15856 error_message <<
"imesh " <<
i_mesh
15857 <<
" is greater than the number of sub meshes "
15869 mmesh_pt->p_refine_uniformly(doc_info);
15873 oomph_info <<
"Info/Warning: Mesh cannot be refined uniformly "
15916 oomph_info <<
"Info/Warning: Mesh cannot be unrefined uniformly "
15976 std::ostringstream error_message;
15977 error_message <<
"imesh " <<
i_mesh
15978 <<
" is greater than the number of sub meshes "
15994 oomph_info <<
"Info/Warning: Mesh cannot be unrefined uniformly "
16037 mmesh_pt->p_unrefine_uniformly(doc_info);
16041 oomph_info <<
"Info/Warning: Mesh cannot be p-unrefined uniformly "
16049 throw OomphLibError(
"This functionality has not yet been tested.",
16059 mmesh_pt->p_unrefine_uniformly(doc_info);
16090 std::ostringstream error_message;
16091 error_message <<
"imesh " <<
i_mesh
16092 <<
" is greater than the number of sub meshes "
16104 mmesh_pt->p_unrefine_uniformly(doc_info);
16108 oomph_info <<
"Info/Warning: Mesh cannot be p-unrefined uniformly "
16148 <<
"\n\n===========================================================\n";
16149 oomph_info <<
" ******** WARNING *********** \n";
16151 <<
"===========================================================\n";
16152 oomph_info <<
"Problem::unsteady_newton_solve() called with "
16156 oomph_info <<
"This doesn't make sense (shifting does have to be done"
16158 oomph_info <<
"since we're constantly re-assigning the initial conditions"
16161 <<
"\n===========================================================\n\n";
16184#ifdef OOMPH_HAS_MPI
16209 <<
n_unrefined <<
" were unrefined, in total." << std::endl;
16214 oomph_info <<
"\n \n Solution is fully converged in "
16215 <<
"Problem::unsteady_newton_solver() \n \n ";
16230 oomph_info <<
"Re-setting initial condition " << std::endl;
16265 <<
"----------------------------------------------------------"
16267 <<
"Reached max. number of adaptations in \n"
16268 <<
"Problem::unsteady_newton_solver().\n"
16269 <<
"----------------------------------------------------------"
16304#ifdef OOMPH_HAS_MPI
16330#ifdef OOMPH_HAS_MPI
16331 <<
", in total (over all processors).\n";
16340 oomph_info <<
"\n \n Solution is fully converged in "
16341 <<
"Problem::newton_solver(). \n \n ";
16363 <<
"USER-DEFINED ERROR IN NEWTON SOLVER " << std::endl;
16368 <<
") REACHED WITHOUT CONVERGENCE " << std::endl;
16380 error_stream <<
"Error occured in adaptive Newton solver. "
16398 <<
"----------------------------------------------------------"
16400 <<
"Reached max. number of adaptations in \n"
16401 <<
"Problem::newton_solver().\n"
16402 <<
"----------------------------------------------------------"
16434#ifdef OOMPH_HAS_MPI
16477 oomph_info <<
"Checking halo schemes on single mesh" << std::endl;
16478 doc_info.
label() =
"_one_and_only_mesh_";
16488 std::stringstream
tmp;
16586 for (
unsigned n = 0;
n <
n_nod;
n++)
16712 for (
unsigned n = 0;
n <
n_nod;
n++)
16823 double t_end = 0.0;
16849 for (
unsigned e = 0;
e <
nelem;
e++)
16874 for (
unsigned j = 0;
j <
nnod;
j++)
16935 oomph_info <<
"Time for copy_haloed_eqn_numbers_helper for halos: "
16949 <<
"Time for copy_haloed_eqn_numbers_helper for external halos: "
16958 if (assign_local_eqn_numbers)
16980 oomph_info <<
"Time for assign_local_eqn_numbers in sync: "
17073 for (
unsigned n = 0;
n <
n_nod;
n++)
17097 ->add_eqn_numbers_to_vector(
send_data);
17200 for (
unsigned n = 0;
n <
n_nod;
n++)
17269 warn_message <<
"WARNING: You've tried to load balance a problem over\n"
17270 <<
"only one processor: ignoring your request.\n";
17272 "Problem::load_balance()",
17284 error_stream <<
"You have called Problem::load_balance()\n"
17285 <<
"on a non-distributed problem. This doesn't\n"
17286 <<
"make sense -- go distribute your problem first."
17394 std::map<GeneralisedElement*, unsigned>
17481 <<
"The number of nonhalo elements (" <<
count_td
17482 <<
") found in (all)\n"
17483 <<
"the (sub)-mesh(es) is different from the number of target "
17486 <<
"Please ensure that you called the rebuild_global_mesh() method "
17487 <<
"after the\npossible deletion of FaceElements in "
17488 <<
"actions_before_distribute()!!!\n\n";
17490 "Problem::load_balance()",
17656 oomph_info <<
"CPU for partition calculation for roots: "
17679 ref_mesh_pt->uniform_refinement_level_when_pruned();
17705 ref_mesh_pt->uniform_refinement_level_when_pruned();
17735 for (
unsigned i = 0;
i <
n;
i++)
17793 for (
unsigned i = 0;
i <
nref;
i++)
17836 <<
nref <<
" times\n";
17837 for (
unsigned i = 0;
i <
nref;
i++)
17877 unsigned count = 0;
17901 <<
"The number of READ target domains for nonhalo elements\n"
17902 <<
" is (" <<
count <<
"), but the number of target domains for\n"
17906 "Problem::load_balance()",
17926 for (
unsigned e = 0;
e <
n_ele;
e++)
17960 for (
unsigned e = 0;
e <
n_ele;
e++)
17976 <<
") is not the same as the number of\nadded elements ("
17977 <<
counter_base <<
") to the Base_mesh_element_pt data "
17978 <<
"structure!!!\n\n";
17980 "Problem::load_balance()",
18007 std::stringstream
info;
18008 info <<
"================================================\n";
18009 info <<
"INFO: I've come across a FaceElement while \n";
18010 info <<
" load-balancing a problem. \n";
18011 info <<
"================================================\n";
18019 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
18043 error_stream <<
"Distributing one-and-only mesh containing "
18054 oomph_info <<
"Distributing one and only mesh\n"
18055 <<
"------------------------------" << std::endl;
18063 new_domain_for_base_element,
18085 <<
n_mesh <<
" in total\n"
18086 <<
"---------------------------------------------"
18097 submesh_element_partition[
i_mesh],
18120 for (
unsigned e = 0;
e <
n_del;
e++)
18161 for (
unsigned e = 0;
e <
n_del;
e++)
18176 oomph_info <<
"CPU for build and distribution of new mesh(es): "
18210 oomph_info <<
"CPU for refinement of base mesh: "
18256 this->
mesh_pt() = old_mesh_pt[0];
18262 error_stream <<
"The only one mesh in the problem is not an "
18263 "unstructured mesh,\n"
18264 <<
"but the flag 'are_there_unstructures_meshes' ("
18266 <<
") was turned on,\n"
18267 <<
"this is weird. Please check for any condition "
18269 <<
"turned on this flag!!!!\n\n";
18271 "Problem::load_balance()",
18335 oomph_info <<
"CPU for transferring solution to new mesh(es): "
18354 <<
"Total number of elements on this processor after load balance: "
18357 oomph_info <<
"Number of non-halo elements on this processor after "
18367 <<
"Number of dofs in load_balance() has changed from " <<
old_ndof
18368 <<
" to " <<
n_dof <<
"\n"
18369 <<
"Check that you've implemented any necessary "
18370 "actions_before/after\n"
18371 <<
"adapt/distribute functions, e.g. to pin redundant pressure dofs"
18445 for (
unsigned h = 0; h <
nhaloed; h++)
18452 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
18477 unsigned count = 0;
18504 for (
unsigned j = 0;
j <
nhalo;
j++)
18561 std::stringstream error_message;
18563 <<
" doesn't match that actually send: "
18716 for (
unsigned j = 0;
j <
nhalo;
j++)
18799 std::stringstream error_message;
18801 <<
" doesn't match that actually send: "
18839 for (
unsigned jd = 0;
jd <
nd;
jd++)
18859 for (
unsigned j = 0;
j <
n;
j++)
18991 for (
unsigned j = 0;
j <
n;
j++)
19019 const int n_proc = comm_pt->nproc();
19092 for (
unsigned b = 0; b <
nbatch; b++)
19125 std::ostringstream error_message;
19127 <<
"Number of leaves: " <<
n_leaf <<
" "
19128 <<
" doesn't match number of elements sent in batch: " <<
nel
19150 std::ostringstream error_message;
19152 <<
"Non-refineable root element should only be associated with"
19153 <<
" one element but nel=" <<
nel <<
"\n";
19163 for (
unsigned e = 0;
e <
nel;
e++)
19174 for (
unsigned j = 0;
j <
nnod;
j++)
19192 std::ostringstream error_message;
19194 <<
"Node has more values, namely " <<
nod_pt->nvalue()
19195 <<
", than we're about to receive, namely " <<
nval
19196 <<
". Something's wrong!\n";
19219 std::ostringstream error_message;
19220 error_message <<
"Local node is boundary node but "
19221 "information sent is\n"
19222 <<
"for non-boundary node\n";
19239 ->index_of_first_value_assigned_by_face_element_pt() ==
19243 ->index_of_first_value_assigned_by_face_element_pt() =
19244 new std::map<unsigned, unsigned>;
19249 std::map<unsigned, unsigned>*
map_pt =
19251 ->index_of_first_value_assigned_by_face_element_pt();
19274 std::ostringstream error_message;
19275 error_message <<
"Local node is not a boundary node but "
19277 <<
"sent is for boundary node.\n";
19353 const int n_proc = comm_pt->nproc();
19374 unsigned max_refinement_level = 0;
19395 std::map<unsigned, Vector<unsigned>>
19444 for (
unsigned e = 0;
e <
nele;
e++)
19476 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
19513 "Base_mesh_element_number_plus_one[...]=0",
19535 root_el_pt->tree_pt()->stick_leaves_into_vector(
19559 if (level > max_refinement_level)
19561 max_refinement_level = level;
19575 std::ostringstream error_message;
19577 <<
"All elements associated with same root must have "
19578 <<
"same target. during load balancing\n";
19595 std::ostringstream error_message;
19599 <<
" target domains for local non-halo elelemts. \n "
19600 <<
"Very Odd -- we do (now) strip out the information for elements\n"
19601 <<
"that are removed in actions_before_distribute()...\n";
19624 comm_pt->mpi_comm());
19640 comm_pt->mpi_comm());
19653 comm_pt->mpi_comm());
19664 std::ostringstream error_message;
19665 error_message <<
"Old domain for base element " <<
j <<
": "
19667 <<
"or its incarnation as refineable el: "
19670 <<
" which is of type "
19673 <<
"appear to have been assigned by any processor\n";
19683 std::ostringstream error_message;
19684 error_message <<
"New domain for base element " <<
j
19685 <<
"which is of type "
19688 <<
"appear to have been assigned by any processor\n";
19729 for (
unsigned b = 0; b <
nbatch; b++)
19745 for (
unsigned e = 0;
e <
nel;
e++)
19756 for (
unsigned j = 0;
j <
nnod;
j++)
19763 unsigned nt =
nod_pt->ntstorage();
19769 for (
unsigned t = 0;
t <
nt;
t++)
19777 for (
unsigned i = 0;
i <
n_dim;
i++)
19816 std::map<unsigned, unsigned>*
map_pt =
19817 bnod_pt->index_of_first_value_assigned_by_face_element_pt();
19831 for (std::map<unsigned, unsigned>::iterator
p =
19836 send_data.push_back(
double((*p).first));
19837 send_data.push_back(
double((*p).second));
19862 std::ostringstream error_message;
19865 <<
" doesn't match total number of elements sent in batch: "
19943 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
19968 "Base_mesh_element_number_plus_one[...]=0",
19977 root_el_pt->tree_pt()->stick_all_tree_nodes_into_vector(
20102 throw OomphLibError(
"Base_mesh_element_number_plus_one[...]=0",
20217 for (
unsigned e = 0;
e <
nel;
e++)
20294 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....
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
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....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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)....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
//////////////////////////////////////////////////////////////////////
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
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...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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...