double diffusion_coefficient;
double reaction_coefficient;
- bool time_dependent_forcing;
+ bool time_dependent_forcing;
unsigned int initial_refinement_level;
unsigned int max_refinement_level;
unsigned int fe_order;
- double start_time;
- double stop_time;
+ double start_time;
+ double stop_time;
unsigned int n_time_steps;
unsigned int save_interval;
unsigned int patch_level;
- void read_parameter_file(const std::string &file_name);
+ void
+ read_parameter_file(const std::string &file_name);
+
private:
- void configure_parameter_handler(ParameterHandler ¶meter_handler);
+ void
+ configure_parameter_handler(ParameterHandler ¶meter_handler);
};
-}
+} // namespace CDR
#endif
{
using namespace dealii;
- template<int dim, typename MatrixType>
- void create_system_matrix
- (const DoFHandler<dim> &dof_handler,
- const QGauss<dim> &quad,
- const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- MatrixType &system_matrix);
+ template <int dim, typename MatrixType>
+ void
+ create_system_matrix(
+ const DoFHandler<dim> & dof_handler,
+ const QGauss<dim> & quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ MatrixType & system_matrix);
- template<int dim, typename MatrixType>
- void create_system_matrix
- (const DoFHandler<dim> &dof_handler,
- const QGauss<dim> &quad,
- const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- const AffineConstraints<double> &constraints,
- MatrixType &system_matrix);
-}
+ template <int dim, typename MatrixType>
+ void
+ create_system_matrix(
+ const DoFHandler<dim> & dof_handler,
+ const QGauss<dim> & quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ const AffineConstraints<double> & constraints,
+ MatrixType & system_matrix);
+} // namespace CDR
#endif
#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/lac/affine_constraints.h>
-
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/lac/affine_constraints.h>
+
#include <deal.II-cdr/parameters.h>
#include <deal.II-cdr/system_matrix.h>
// This is the actual implementation of the <code>create_system_matrix</code>
// function described in the header file. It is similar to the system matrix
// assembly routine in step-40.
- template<int dim, typename UpdateFunction>
- void internal_create_system_matrix
- (const DoFHandler<dim> &dof_handler,
- const QGauss<dim> &quad,
- const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- UpdateFunction update_system_matrix)
+ template <int dim, typename UpdateFunction>
+ void
+ internal_create_system_matrix(
+ const DoFHandler<dim> & dof_handler,
+ const QGauss<dim> & quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ UpdateFunction update_system_matrix)
{
- auto &fe = dof_handler.get_fe();
- const auto dofs_per_cell = fe.dofs_per_cell;
+ auto & fe = dof_handler.get_fe();
+ const auto dofs_per_cell = fe.dofs_per_cell;
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
- FEValues<dim> fe_values(fe, quad, update_values | update_gradients |
- update_quadrature_points | update_JxW_values);
+ FEValues<dim> fe_values(fe,
+ quad,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
std::vector<types::global_dof_index> local_indices(dofs_per_cell);
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
- const auto convection_contribution = current_convection
- *fe_values.shape_grad(j, q);
- cell_matrix(i, j) += fe_values.JxW(q)*
- // Here are the time step, mass, and reaction parts:
- ((1.0 + time_step/2.0*parameters.reaction_coefficient)
- *fe_values.shape_value(i, q)*fe_values.shape_value(j, q)
- + time_step/2.0*
- // and the convection part:
- (fe_values.shape_value(i, q)*convection_contribution
- // and, finally, the diffusion part:
- + parameters.diffusion_coefficient
- *(fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q)))
- );
+ const auto convection_contribution =
+ current_convection * fe_values.shape_grad(j, q);
+ cell_matrix(i, j) +=
+ fe_values.JxW(q) *
+ // Here are the time step, mass, and reaction parts:
+ ((1.0 +
+ time_step / 2.0 * parameters.reaction_coefficient) *
+ fe_values.shape_value(i, q) *
+ fe_values.shape_value(j, q) +
+ time_step / 2.0 *
+ // and the convection part:
+ (fe_values.shape_value(i, q) *
+ convection_contribution
+ // and, finally, the diffusion part:
+ + parameters.diffusion_coefficient *
+ (fe_values.shape_grad(i, q) *
+ fe_values.shape_grad(j, q))));
}
}
}
}
}
- template<int dim, typename MatrixType>
- void create_system_matrix
- (const DoFHandler<dim> &dof_handler,
- const QGauss<dim> &quad,
- const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- const AffineConstraints<double> &constraints,
- MatrixType &system_matrix)
+ template <int dim, typename MatrixType>
+ void
+ create_system_matrix(
+ const DoFHandler<dim> & dof_handler,
+ const QGauss<dim> & quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ const AffineConstraints<double> & constraints,
+ MatrixType & system_matrix)
{
- internal_create_system_matrix<dim>
- (dof_handler, quad, convection_function, parameters, time_step,
- [&constraints, &system_matrix](const std::vector<types::global_dof_index> &local_indices,
- const FullMatrix<double> &cell_matrix)
- {
- constraints.distribute_local_to_global
- (cell_matrix, local_indices, system_matrix);
- });
+ internal_create_system_matrix<dim>(
+ dof_handler,
+ quad,
+ convection_function,
+ parameters,
+ time_step,
+ [&constraints, &system_matrix](
+ const std::vector<types::global_dof_index> &local_indices,
+ const FullMatrix<double> & cell_matrix) {
+ constraints.distribute_local_to_global(cell_matrix,
+ local_indices,
+ system_matrix);
+ });
}
- template<int dim, typename MatrixType>
- void create_system_matrix
- (const DoFHandler<dim> &dof_handler,
- const QGauss<dim> &quad,
- const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- MatrixType &system_matrix)
+ template <int dim, typename MatrixType>
+ void
+ create_system_matrix(
+ const DoFHandler<dim> & dof_handler,
+ const QGauss<dim> & quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ MatrixType & system_matrix)
{
- internal_create_system_matrix<dim>
- (dof_handler, quad, convection_function, parameters, time_step,
- [&system_matrix](const std::vector<types::global_dof_index> &local_indices,
- const FullMatrix<double> &cell_matrix)
- {
- system_matrix.add(local_indices, cell_matrix);
- });
+ internal_create_system_matrix<dim>(
+ dof_handler,
+ quad,
+ convection_function,
+ parameters,
+ time_step,
+ [&system_matrix](
+ const std::vector<types::global_dof_index> &local_indices,
+ const FullMatrix<double> & cell_matrix) {
+ system_matrix.add(local_indices, cell_matrix);
+ });
}
-}
+} // namespace CDR
#endif
{
using namespace dealii;
- template<int dim, typename VectorType>
- void create_system_rhs
- (const DoFHandler<dim> &dof_handler,
- const QGauss<dim> &quad,
- const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
- const std::function<double(double, const Point<dim>)> &forcing_function,
- const CDR::Parameters ¶meters,
- const VectorType &previous_solution,
- const AffineConstraints<double> &constraints,
- const double current_time,
- VectorType &system_rhs);
-}
+ template <int dim, typename VectorType>
+ void
+ create_system_rhs(
+ const DoFHandler<dim> & dof_handler,
+ const QGauss<dim> & quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+ const std::function<double(double, const Point<dim>)> &forcing_function,
+ const CDR::Parameters & parameters,
+ const VectorType & previous_solution,
+ const AffineConstraints<double> & constraints,
+ const double current_time,
+ VectorType & system_rhs);
+} // namespace CDR
#endif
{
using namespace dealii;
- template<int dim, typename VectorType>
- void create_system_rhs
- (const DoFHandler<dim> &dof_handler,
- const QGauss<dim> &quad,
- const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
- const std::function<double(double, const Point<dim>)> &forcing_function,
- const CDR::Parameters ¶meters,
- const VectorType &previous_solution,
- const AffineConstraints<double> &constraints,
- const double current_time,
- VectorType &system_rhs)
+ template <int dim, typename VectorType>
+ void
+ create_system_rhs(
+ const DoFHandler<dim> & dof_handler,
+ const QGauss<dim> & quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+ const std::function<double(double, const Point<dim>)> &forcing_function,
+ const CDR::Parameters & parameters,
+ const VectorType & previous_solution,
+ const AffineConstraints<double> & constraints,
+ const double current_time,
+ VectorType & system_rhs)
{
- auto &fe = dof_handler.get_fe();
- const auto dofs_per_cell = fe.dofs_per_cell;
- const double time_step = (parameters.stop_time - parameters.start_time)
- /parameters.n_time_steps;
- FEValues<dim> fe_values(fe, quad, update_values | update_gradients |
- update_quadrature_points | update_JxW_values);
+ auto & fe = dof_handler.get_fe();
+ const auto dofs_per_cell = fe.dofs_per_cell;
+ const double time_step =
+ (parameters.stop_time - parameters.start_time) / parameters.n_time_steps;
+ FEValues<dim> fe_values(fe,
+ quad,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
- Vector<double> cell_rhs(dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
- Vector<double> current_fe_coefficients(dofs_per_cell);
+ Vector<double> current_fe_coefficients(dofs_per_cell);
std::vector<types::global_dof_index> local_indices(dofs_per_cell);
- const double previous_time {current_time - time_step};
+ const double previous_time{current_time - time_step};
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell->get_dof_indices(local_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
- current_fe_coefficients[i] = previous_solution[local_indices[i]];
+ current_fe_coefficients[i] =
+ previous_solution[local_indices[i]];
}
for (unsigned int q = 0; q < quad.size(); ++q)
const auto current_convection =
convection_function(fe_values.quadrature_point(q));
- const double current_forcing = forcing_function
- (current_time, fe_values.quadrature_point(q));
- const double previous_forcing = forcing_function
- (previous_time, fe_values.quadrature_point(q));
+ const double current_forcing =
+ forcing_function(current_time, fe_values.quadrature_point(q));
+ const double previous_forcing =
+ forcing_function(previous_time,
+ fe_values.quadrature_point(q));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
- const auto convection_contribution = current_convection
- *fe_values.shape_grad(j, q);
+ const auto convection_contribution =
+ current_convection * fe_values.shape_grad(j, q);
- cell_rhs(i) += fe_values.JxW(q)*
- // Here are the mass and reaction part:
- (((1.0 - time_step/2.0*parameters.reaction_coefficient)
- *fe_values.shape_value(i, q)*fe_values.shape_value(j, q)
- - time_step/2.0*
- // the convection part:
- (fe_values.shape_value(i, q)*convection_contribution
- // the diffusion part:
- + parameters.diffusion_coefficient
- *(fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q))))
- *current_fe_coefficients[j]
- // and, finally, the forcing function part:
- + time_step/2.0*
- (current_forcing + previous_forcing)
- *fe_values.shape_value(i, q));
+ cell_rhs(i) +=
+ fe_values.JxW(q) *
+ // Here are the mass and reaction part:
+ (((1.0 - time_step / 2.0 *
+ parameters.reaction_coefficient) *
+ fe_values.shape_value(i, q) *
+ fe_values.shape_value(j, q) -
+ time_step / 2.0 *
+ // the convection part:
+ (fe_values.shape_value(i, q) *
+ convection_contribution
+ // the diffusion part:
+ + parameters.diffusion_coefficient *
+ (fe_values.shape_grad(i, q) *
+ fe_values.shape_grad(j, q)))) *
+ current_fe_coefficients[j]
+ // and, finally, the forcing function part:
+ + time_step / 2.0 *
+ (current_forcing + previous_forcing) *
+ fe_values.shape_value(i, q));
}
}
}
- constraints.distribute_local_to_global(cell_rhs, local_indices, system_rhs);
+ constraints.distribute_local_to_global(cell_rhs,
+ local_indices,
+ system_rhs);
}
}
}
-}
+} // namespace CDR
#endif
public:
WritePVTUOutput(const unsigned int patch_level);
- template<int dim, typename VectorType>
- void write_output(const DoFHandler<dim> &dof_handler,
- const VectorType &solution,
- const unsigned int time_step_n,
- const double current_time);
+ template <int dim, typename VectorType>
+ void
+ write_output(const DoFHandler<dim> &dof_handler,
+ const VectorType & solution,
+ const unsigned int time_step_n,
+ const double current_time);
+
private:
const unsigned int patch_level;
const unsigned int this_mpi_process;
};
-}
+} // namespace CDR
#endif
#include <deal.II-cdr/write_pvtu_output.h>
-#include <string>
#include <fstream>
+#include <string>
#include <vector>
// Here is the implementation of the important function. This is similar to
{
using namespace dealii;
- template<int dim, typename VectorType>
- void WritePVTUOutput::write_output(const DoFHandler<dim> &dof_handler,
- const VectorType &solution,
- const unsigned int time_step_n,
- const double current_time)
+ template <int dim, typename VectorType>
+ void
+ WritePVTUOutput::write_output(const DoFHandler<dim> &dof_handler,
+ const VectorType & solution,
+ const unsigned int time_step_n,
+ const double current_time)
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "u");
- const auto &triangulation = dof_handler.get_triangulation();
- Vector<float> subdomain (triangulation.n_active_cells());
+ const auto & triangulation = dof_handler.get_triangulation();
+ Vector<float> subdomain(triangulation.n_active_cells());
for (auto &domain : subdomain)
{
domain = triangulation.locally_owned_subdomain();
flags.time = current_time;
// While the default flag is for the best compression level, using
// <code>best_speed</code> makes this function much faster.
- flags.compression_level = DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
+ flags.compression_level =
+ DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
data_out.set_flags(flags);
unsigned int subdomain_n;
subdomain_n = triangulation.locally_owned_subdomain();
}
- std::ofstream output
- ("solution-" + Utilities::int_to_string(time_step_n) + "."
- + Utilities::int_to_string(subdomain_n, 4)
- + ".vtu");
+ std::ofstream output("solution-" + Utilities::int_to_string(time_step_n) +
+ "." + Utilities::int_to_string(subdomain_n, 4) +
+ ".vtu");
data_out.write_vtu(output);
if (this_mpi_process == 0)
{
std::vector<std::string> filenames;
- for (unsigned int i = 0; i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+ for (unsigned int i = 0;
+ i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
++i)
- filenames.push_back
- ("solution-" + Utilities::int_to_string (time_step_n) + "."
- + Utilities::int_to_string (i, 4) + ".vtu");
- std::ofstream master_output
- ("solution-" + Utilities::int_to_string(time_step_n) + ".pvtu");
+ filenames.push_back("solution-" +
+ Utilities::int_to_string(time_step_n) + "." +
+ Utilities::int_to_string(i, 4) + ".vtu");
+ std::ofstream master_output(
+ "solution-" + Utilities::int_to_string(time_step_n) + ".pvtu");
data_out.write_pvtu_record(master_output, filenames);
}
}
-}
+} // namespace CDR
#endif
namespace CDR
{
- void Parameters::configure_parameter_handler(ParameterHandler ¶meter_handler)
+ void
+ Parameters::configure_parameter_handler(ParameterHandler ¶meter_handler)
{
parameter_handler.enter_subsection("Geometry");
{
- parameter_handler.declare_entry
- ("inner_radius", "1.0", Patterns::Double(0.0), "Inner radius.");
- parameter_handler.declare_entry
- ("outer_radius", "2.0", Patterns::Double(0.0), "Outer radius.");
+ parameter_handler.declare_entry("inner_radius",
+ "1.0",
+ Patterns::Double(0.0),
+ "Inner radius.");
+ parameter_handler.declare_entry("outer_radius",
+ "2.0",
+ Patterns::Double(0.0),
+ "Outer radius.");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Physical Parameters");
{
- parameter_handler.declare_entry
- ("diffusion_coefficient", "1.0", Patterns::Double(0.0), "Diffusion coefficient.");
- parameter_handler.declare_entry
- ("reaction_coefficient", "1.0", Patterns::Double(0.0), "Reaction coefficient.");
- parameter_handler.declare_entry
- ("time_dependent_forcing", "true", Patterns::Bool(), "Whether or not "
- "the forcing function depends on time.");
+ parameter_handler.declare_entry("diffusion_coefficient",
+ "1.0",
+ Patterns::Double(0.0),
+ "Diffusion coefficient.");
+ parameter_handler.declare_entry("reaction_coefficient",
+ "1.0",
+ Patterns::Double(0.0),
+ "Reaction coefficient.");
+ parameter_handler.declare_entry("time_dependent_forcing",
+ "true",
+ Patterns::Bool(),
+ "Whether or not "
+ "the forcing function depends on time.");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Finite Element");
{
- parameter_handler.declare_entry
- ("initial_refinement_level", "1", Patterns::Integer(1),
- "Initial number of levels in the mesh.");
- parameter_handler.declare_entry
- ("max_refinement_level", "1", Patterns::Integer(1),
- "Maximum number of levels in the mesh.");
- parameter_handler.declare_entry
- ("fe_order", "1", Patterns::Integer(1), "Finite element order.");
+ parameter_handler.declare_entry("initial_refinement_level",
+ "1",
+ Patterns::Integer(1),
+ "Initial number of levels in the mesh.");
+ parameter_handler.declare_entry("max_refinement_level",
+ "1",
+ Patterns::Integer(1),
+ "Maximum number of levels in the mesh.");
+ parameter_handler.declare_entry("fe_order",
+ "1",
+ Patterns::Integer(1),
+ "Finite element order.");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Time Step");
{
- parameter_handler.declare_entry
- ("start_time", "0.0", Patterns::Double(0.0), "Start time.");
- parameter_handler.declare_entry
- ("stop_time", "1.0", Patterns::Double(1.0), "Stop time.");
- parameter_handler.declare_entry
- ("n_time_steps", "1", Patterns::Integer(1), "Number of time steps.");
+ parameter_handler.declare_entry("start_time",
+ "0.0",
+ Patterns::Double(0.0),
+ "Start time.");
+ parameter_handler.declare_entry("stop_time",
+ "1.0",
+ Patterns::Double(1.0),
+ "Stop time.");
+ parameter_handler.declare_entry("n_time_steps",
+ "1",
+ Patterns::Integer(1),
+ "Number of time steps.");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Output");
{
- parameter_handler.declare_entry
- ("save_interval", "10", Patterns::Integer(1), "Save interval.");
- parameter_handler.declare_entry
- ("patch_level", "2", Patterns::Integer(0), "Patch level.");
+ parameter_handler.declare_entry("save_interval",
+ "10",
+ Patterns::Integer(1),
+ "Save interval.");
+ parameter_handler.declare_entry("patch_level",
+ "2",
+ Patterns::Integer(0),
+ "Patch level.");
}
parameter_handler.leave_subsection();
}
- void Parameters::read_parameter_file(const std::string &file_name)
+ void
+ Parameters::read_parameter_file(const std::string &file_name)
{
ParameterHandler parameter_handler;
{
parameter_handler.enter_subsection("Physical Parameters");
{
- diffusion_coefficient = parameter_handler.get_double("diffusion_coefficient");
- reaction_coefficient = parameter_handler.get_double("reaction_coefficient");
- time_dependent_forcing = parameter_handler.get_bool("time_dependent_forcing");
+ diffusion_coefficient =
+ parameter_handler.get_double("diffusion_coefficient");
+ reaction_coefficient =
+ parameter_handler.get_double("reaction_coefficient");
+ time_dependent_forcing =
+ parameter_handler.get_bool("time_dependent_forcing");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Finite Element");
{
- initial_refinement_level = parameter_handler.get_integer("initial_refinement_level");
- max_refinement_level = parameter_handler.get_integer("max_refinement_level");
+ initial_refinement_level =
+ parameter_handler.get_integer("initial_refinement_level");
+ max_refinement_level =
+ parameter_handler.get_integer("max_refinement_level");
fe_order = parameter_handler.get_integer("fe_order");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Time Step");
{
- start_time = parameter_handler.get_double("start_time");
- stop_time = parameter_handler.get_double("stop_time");
+ start_time = parameter_handler.get_double("start_time");
+ stop_time = parameter_handler.get_double("stop_time");
n_time_steps = parameter_handler.get_integer("n_time_steps");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Output");
{
save_interval = parameter_handler.get_integer("save_interval");
- patch_level = parameter_handler.get_integer("patch_level");
+ patch_level = parameter_handler.get_integer("patch_level");
}
parameter_handler.leave_subsection();
}
-}
+} // namespace CDR
{
using namespace dealii;
- template
- void create_system_matrix<2, SparseMatrix<double>>
- (const DoFHandler<2> &dof_handler,
- const QGauss<2> &quad,
- const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- SparseMatrix<double> &system_matrix);
+ template void
+ create_system_matrix<2, SparseMatrix<double>>(
+ const DoFHandler<2> & dof_handler,
+ const QGauss<2> & quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ SparseMatrix<double> & system_matrix);
- template
- void create_system_matrix<3, SparseMatrix<double>>
- (const DoFHandler<3> &dof_handler,
- const QGauss<3> &quad,
- const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- SparseMatrix<double> &system_matrix);
+ template void
+ create_system_matrix<3, SparseMatrix<double>>(
+ const DoFHandler<3> & dof_handler,
+ const QGauss<3> & quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ SparseMatrix<double> & system_matrix);
- template
- void create_system_matrix<2, SparseMatrix<double>>
- (const DoFHandler<2> &dof_handler,
- const QGauss<2> &quad,
- const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- const AffineConstraints<double> &constraints,
- SparseMatrix<double> &system_matrix);
+ template void
+ create_system_matrix<2, SparseMatrix<double>>(
+ const DoFHandler<2> & dof_handler,
+ const QGauss<2> & quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ const AffineConstraints<double> & constraints,
+ SparseMatrix<double> & system_matrix);
- template
- void create_system_matrix<3, SparseMatrix<double>>
- (const DoFHandler<3> &dof_handler,
- const QGauss<3> &quad,
- const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- const AffineConstraints<double> &constraints,
- SparseMatrix<double> &system_matrix);
+ template void
+ create_system_matrix<3, SparseMatrix<double>>(
+ const DoFHandler<3> & dof_handler,
+ const QGauss<3> & quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ const AffineConstraints<double> & constraints,
+ SparseMatrix<double> & system_matrix);
- template
- void create_system_matrix<2, TrilinosWrappers::SparseMatrix>
- (const DoFHandler<2> &dof_handler,
- const QGauss<2> &quad,
- const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- TrilinosWrappers::SparseMatrix &system_matrix);
+ template void
+ create_system_matrix<2, TrilinosWrappers::SparseMatrix>(
+ const DoFHandler<2> & dof_handler,
+ const QGauss<2> & quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ TrilinosWrappers::SparseMatrix & system_matrix);
- template
- void create_system_matrix<3, TrilinosWrappers::SparseMatrix>
- (const DoFHandler<3> &dof_handler,
- const QGauss<3> &quad,
- const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- TrilinosWrappers::SparseMatrix &system_matrix);
+ template void
+ create_system_matrix<3, TrilinosWrappers::SparseMatrix>(
+ const DoFHandler<3> & dof_handler,
+ const QGauss<3> & quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ TrilinosWrappers::SparseMatrix & system_matrix);
- template
- void create_system_matrix<2, TrilinosWrappers::SparseMatrix>
- (const DoFHandler<2> &dof_handler,
- const QGauss<2> &quad,
- const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- const AffineConstraints<double> &constraints,
- TrilinosWrappers::SparseMatrix &system_matrix);
+ template void
+ create_system_matrix<2, TrilinosWrappers::SparseMatrix>(
+ const DoFHandler<2> & dof_handler,
+ const QGauss<2> & quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ const AffineConstraints<double> & constraints,
+ TrilinosWrappers::SparseMatrix & system_matrix);
- template
- void create_system_matrix<3, TrilinosWrappers::SparseMatrix>
- (const DoFHandler<3> &dof_handler,
- const QGauss<3> &quad,
- const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
- const CDR::Parameters ¶meters,
- const double time_step,
- const AffineConstraints<double> &constraints,
- TrilinosWrappers::SparseMatrix &system_matrix);
-}
+ template void
+ create_system_matrix<3, TrilinosWrappers::SparseMatrix>(
+ const DoFHandler<3> & dof_handler,
+ const QGauss<3> & quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
+ const CDR::Parameters & parameters,
+ const double time_step,
+ const AffineConstraints<double> & constraints,
+ TrilinosWrappers::SparseMatrix & system_matrix);
+} // namespace CDR
{
using namespace dealii;
- template
- void create_system_rhs<2, Vector<double>>
- (const DoFHandler<2> &dof_handler,
- const QGauss<2> &quad,
- const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
- const std::function<double(double, const Point<2>)> &forcing_function,
- const CDR::Parameters ¶meters,
- const Vector<double> &previous_solution,
- const AffineConstraints<double> &constraints,
- const double current_time,
- Vector<double> &system_rhs);
+ template void
+ create_system_rhs<2, Vector<double>>(
+ const DoFHandler<2> & dof_handler,
+ const QGauss<2> & quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> & convection_function,
+ const std::function<double(double, const Point<2>)> &forcing_function,
+ const CDR::Parameters & parameters,
+ const Vector<double> & previous_solution,
+ const AffineConstraints<double> & constraints,
+ const double current_time,
+ Vector<double> & system_rhs);
- template
- void create_system_rhs<3, Vector<double>>
- (const DoFHandler<3> &dof_handler,
- const QGauss<3> &quad,
- const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
- const std::function<double(double, const Point<3>)> &forcing_function,
- const CDR::Parameters ¶meters,
- const Vector<double> &previous_solution,
- const AffineConstraints<double> &constraints,
- const double current_time,
- Vector<double> &system_rhs);
+ template void
+ create_system_rhs<3, Vector<double>>(
+ const DoFHandler<3> & dof_handler,
+ const QGauss<3> & quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> & convection_function,
+ const std::function<double(double, const Point<3>)> &forcing_function,
+ const CDR::Parameters & parameters,
+ const Vector<double> & previous_solution,
+ const AffineConstraints<double> & constraints,
+ const double current_time,
+ Vector<double> & system_rhs);
- template
- void create_system_rhs<2, TrilinosWrappers::MPI::Vector>
- (const DoFHandler<2> &dof_handler,
- const QGauss<2> &quad,
- const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
- const std::function<double(double, const Point<2>)> &forcing_function,
- const CDR::Parameters ¶meters,
- const TrilinosWrappers::MPI::Vector &previous_solution,
- const AffineConstraints<double> &constraints,
- const double current_time,
- TrilinosWrappers::MPI::Vector &system_rhs);
+ template void
+ create_system_rhs<2, TrilinosWrappers::MPI::Vector>(
+ const DoFHandler<2> & dof_handler,
+ const QGauss<2> & quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> & convection_function,
+ const std::function<double(double, const Point<2>)> &forcing_function,
+ const CDR::Parameters & parameters,
+ const TrilinosWrappers::MPI::Vector & previous_solution,
+ const AffineConstraints<double> & constraints,
+ const double current_time,
+ TrilinosWrappers::MPI::Vector & system_rhs);
- template
- void create_system_rhs<3, TrilinosWrappers::MPI::Vector>
- (const DoFHandler<3> &dof_handler,
- const QGauss<3> &quad,
- const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
- const std::function<double(double, const Point<3>)> &forcing_function,
- const CDR::Parameters ¶meters,
- const TrilinosWrappers::MPI::Vector &previous_solution,
- const AffineConstraints<double> &constraints,
- const double current_time,
- TrilinosWrappers::MPI::Vector &system_rhs);
-}
+ template void
+ create_system_rhs<3, TrilinosWrappers::MPI::Vector>(
+ const DoFHandler<3> & dof_handler,
+ const QGauss<3> & quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> & convection_function,
+ const std::function<double(double, const Point<3>)> &forcing_function,
+ const CDR::Parameters & parameters,
+ const TrilinosWrappers::MPI::Vector & previous_solution,
+ const AffineConstraints<double> & constraints,
+ const double current_time,
+ TrilinosWrappers::MPI::Vector & system_rhs);
+} // namespace CDR
using namespace dealii;
WritePVTUOutput::WritePVTUOutput(const unsigned int patch_level)
- : patch_level {patch_level},
- this_mpi_process {Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)}
+ : patch_level{patch_level}
+ , this_mpi_process{Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)}
{}
- template
- void WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler,
- const Vector<double> &solution,
- const unsigned int time_step_n,
- const double current_time);
+ template void
+ WritePVTUOutput::write_output(const DoFHandler<2> & dof_handler,
+ const Vector<double> &solution,
+ const unsigned int time_step_n,
+ const double current_time);
- template
- void WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler,
- const Vector<double> &solution,
- const unsigned int time_step_n,
- const double current_time);
+ template void
+ WritePVTUOutput::write_output(const DoFHandler<3> & dof_handler,
+ const Vector<double> &solution,
+ const unsigned int time_step_n,
+ const double current_time);
- template
- void WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler,
- const TrilinosWrappers::MPI::Vector &solution,
- const unsigned int time_step_n,
- const double current_time);
+ template void
+ WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler,
+ const TrilinosWrappers::MPI::Vector &solution,
+ const unsigned int time_step_n,
+ const double current_time);
- template
- void WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler,
- const TrilinosWrappers::MPI::Vector &solution,
- const unsigned int time_step_n,
- const double current_time);
-}
+ template void
+ WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler,
+ const TrilinosWrappers::MPI::Vector &solution,
+ const unsigned int time_step_n,
+ const double current_time);
+} // namespace CDR
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/numerics/error_estimator.h>
// These headers are for distributed computations:
-#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
-#include <deal.II/distributed/tria.h>
+#include <deal.II/base/utilities.h>
+
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
-#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
-#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_vector.h>
-#include <chrono>
-#include <functional>
-#include <iostream>
-
+#include <deal.II-cdr/parameters.h>
#include <deal.II-cdr/system_matrix.h>
#include <deal.II-cdr/system_rhs.h>
-#include <deal.II-cdr/parameters.h>
#include <deal.II-cdr/write_pvtu_output.h>
+#include <chrono>
+#include <functional>
+#include <iostream>
+
using namespace dealii;
-constexpr int manifold_id {0};
+constexpr int manifold_id{0};
// This is the actual solver class which performs time iteration and calls the
// appropriate library functions to do it.
-template<int dim>
+template <int dim>
class CDRProblem
{
public:
CDRProblem(const CDR::Parameters ¶meters);
- void run();
+ void
+ run();
+
private:
const CDR::Parameters parameters;
- const double time_step;
- double current_time;
+ const double time_step;
+ double current_time;
- MPI_Comm mpi_communicator;
+ MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
- FE_Q<dim> fe;
- QGauss<dim> quad;
- const SphericalManifold<dim> boundary_description;
+ FE_Q<dim> fe;
+ QGauss<dim> quad;
+ const SphericalManifold<dim> boundary_description;
parallel::distributed::Triangulation<dim> triangulation;
- DoFHandler<dim> dof_handler;
+ DoFHandler<dim> dof_handler;
const std::function<Tensor<1, dim>(const Point<dim>)> convection_function;
const std::function<double(double, const Point<dim>)> forcing_function;
IndexSet locally_relevant_dofs;
AffineConstraints<double> constraints;
- bool first_run;
+ bool first_run;
// As is usual in parallel programs, I keep two copies of parts of the
// complete solution: <code>locally_relevant_solution</code> contains both
// boundary (the @ref GlossGhostCell "ghost cells") while
// <code>completely_distributed_solution</code> only contains the parts of
// the solution computed on the current @ref GlossMPIProcess "MPI process".
- TrilinosWrappers::MPI::Vector locally_relevant_solution;
- TrilinosWrappers::MPI::Vector completely_distributed_solution;
- TrilinosWrappers::MPI::Vector system_rhs;
- TrilinosWrappers::SparseMatrix system_matrix;
+ TrilinosWrappers::MPI::Vector locally_relevant_solution;
+ TrilinosWrappers::MPI::Vector completely_distributed_solution;
+ TrilinosWrappers::MPI::Vector system_rhs;
+ TrilinosWrappers::SparseMatrix system_matrix;
TrilinosWrappers::PreconditionAMG preconditioner;
ConditionalOStream pcout;
- void setup_geometry();
- void setup_system();
- void setup_dofs();
- void refine_mesh();
- void time_iterate();
+ void
+ setup_geometry();
+ void
+ setup_system();
+ void
+ setup_dofs();
+ void
+ refine_mesh();
+ void
+ time_iterate();
};
-template<int dim>
-CDRProblem<dim>::CDRProblem(const CDR::Parameters ¶meters) :
- parameters(parameters),
- time_step {(parameters.stop_time - parameters.start_time)
- /parameters.n_time_steps
-},
-current_time {parameters.start_time},
-mpi_communicator (MPI_COMM_WORLD),
-n_mpi_processes {Utilities::MPI::n_mpi_processes(mpi_communicator)},
-this_mpi_process {Utilities::MPI::this_mpi_process(mpi_communicator)},
-fe(parameters.fe_order),
-quad(parameters.fe_order + 2),
-boundary_description(Point<dim>()),
-triangulation(mpi_communicator, typename Triangulation<dim>::MeshSmoothing
- (Triangulation<dim>::smoothing_on_refinement |
- Triangulation<dim>::smoothing_on_coarsening)),
-dof_handler(triangulation),
-convection_function
-{
- [](const Point<dim> p) -> Tensor<1, dim>
- {Tensor<1, dim> v; v[0] = -p[1]; v[1] = p[0]; return v;}
-},
-forcing_function
-{
- [](double t, const Point<dim> p) -> double
- {
- return std::exp(-8*t)*std::exp(-40*Utilities::fixed_power<6>(p[0] - 1.5))
- *std::exp(-40*Utilities::fixed_power<6>(p[1]));
- }
-},
-first_run {true},
-pcout (std::cout, this_mpi_process == 0)
+template <int dim>
+CDRProblem<dim>::CDRProblem(const CDR::Parameters ¶meters)
+ : parameters(parameters)
+ , time_step{(parameters.stop_time - parameters.start_time) /
+ parameters.n_time_steps}
+ , current_time{parameters.start_time}
+ , mpi_communicator(MPI_COMM_WORLD)
+ , n_mpi_processes{Utilities::MPI::n_mpi_processes(mpi_communicator)}
+ , this_mpi_process{Utilities::MPI::this_mpi_process(mpi_communicator)}
+ , fe(parameters.fe_order)
+ , quad(parameters.fe_order + 2)
+ , boundary_description(Point<dim>())
+ , triangulation(mpi_communicator,
+ typename Triangulation<dim>::MeshSmoothing(
+ Triangulation<dim>::smoothing_on_refinement |
+ Triangulation<dim>::smoothing_on_coarsening))
+ , dof_handler(triangulation)
+ , convection_function{[](const Point<dim> p) -> Tensor<1, dim> {
+ Tensor<1, dim> v;
+ v[0] = -p[1];
+ v[1] = p[0];
+ return v;
+ }}
+ , forcing_function{[](double t, const Point<dim> p) -> double {
+ return std::exp(-8 * t) *
+ std::exp(-40 * Utilities::fixed_power<6>(p[0] - 1.5)) *
+ std::exp(-40 * Utilities::fixed_power<6>(p[1]));
+ }}
+ , first_run{true}
+ , pcout(std::cout, this_mpi_process == 0)
{
Assert(dim == 2, ExcNotImplemented());
}
-template<int dim>
-void CDRProblem<dim>::setup_geometry()
+template <int dim>
+void
+CDRProblem<dim>::setup_geometry()
{
const Point<dim> center;
- GridGenerator::hyper_shell(triangulation, center, parameters.inner_radius,
+ GridGenerator::hyper_shell(triangulation,
+ center,
+ parameters.inner_radius,
parameters.outer_radius);
triangulation.set_manifold(manifold_id, boundary_description);
for (const auto &cell : triangulation.active_cell_iterators())
}
-template<int dim>
-void CDRProblem<dim>::setup_dofs()
+template <int dim>
+void
+CDRProblem<dim>::setup_dofs()
{
dof_handler.distribute_dofs(fe);
- pcout << "Number of degrees of freedom: "
- << dof_handler.n_dofs()
+ pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
- DoFTools::make_zero_boundary_constraints(dof_handler, manifold_id, constraints);
+ DoFTools::make_zero_boundary_constraints(dof_handler,
+ manifold_id,
+ constraints);
constraints.close();
- completely_distributed_solution.reinit
- (locally_owned_dofs, mpi_communicator);
+ completely_distributed_solution.reinit(locally_owned_dofs, mpi_communicator);
- locally_relevant_solution.reinit(locally_owned_dofs, locally_relevant_dofs,
+ locally_relevant_solution.reinit(locally_owned_dofs,
+ locally_relevant_dofs,
mpi_communicator);
}
-template<int dim>
-void CDRProblem<dim>::setup_system()
+template <int dim>
+void
+CDRProblem<dim>::setup_system()
{
DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern(dof_handler, dynamic_sparsity_pattern,
- constraints, /*keep_constrained_dofs*/true);
- SparsityTools::distribute_sparsity_pattern
- (dynamic_sparsity_pattern, dof_handler.locally_owned_dofs(),
- mpi_communicator, locally_relevant_dofs);
+ DoFTools::make_sparsity_pattern(dof_handler,
+ dynamic_sparsity_pattern,
+ constraints,
+ /*keep_constrained_dofs*/ true);
+ SparsityTools::distribute_sparsity_pattern(dynamic_sparsity_pattern,
+ dof_handler.locally_owned_dofs(),
+ mpi_communicator,
+ locally_relevant_dofs);
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
- system_matrix.reinit(locally_owned_dofs, dynamic_sparsity_pattern,
+ system_matrix.reinit(locally_owned_dofs,
+ dynamic_sparsity_pattern,
mpi_communicator);
- CDR::create_system_matrix<dim>
- (dof_handler, quad, convection_function, parameters, time_step, constraints,
- system_matrix);
+ CDR::create_system_matrix<dim>(dof_handler,
+ quad,
+ convection_function,
+ parameters,
+ time_step,
+ constraints,
+ system_matrix);
system_matrix.compress(VectorOperation::add);
preconditioner.initialize(system_matrix);
}
-template<int dim>
-void CDRProblem<dim>::time_iterate()
+template <int dim>
+void
+CDRProblem<dim>::time_iterate()
{
- double current_time = parameters.start_time;
+ double current_time = parameters.start_time;
CDR::WritePVTUOutput pvtu_output(parameters.patch_level);
for (unsigned int time_step_n = 0; time_step_n < parameters.n_time_steps;
++time_step_n)
current_time += time_step;
system_rhs = 0.0;
- CDR::create_system_rhs<dim>
- (dof_handler, quad, convection_function, forcing_function, parameters,
- locally_relevant_solution, constraints, current_time, system_rhs);
+ CDR::create_system_rhs<dim>(dof_handler,
+ quad,
+ convection_function,
+ forcing_function,
+ parameters,
+ locally_relevant_solution,
+ constraints,
+ current_time,
+ system_rhs);
system_rhs.compress(VectorOperation::add);
- SolverControl solver_control(dof_handler.n_dofs(),
- 1e-6*system_rhs.l2_norm(),
+ SolverControl solver_control(dof_handler.n_dofs(),
+ 1e-6 * system_rhs.l2_norm(),
/*log_history = */ false,
/*log_result = */ false);
TrilinosWrappers::SolverGMRES solver(solver_control);
- solver.solve(system_matrix, completely_distributed_solution, system_rhs,
+ solver.solve(system_matrix,
+ completely_distributed_solution,
+ system_rhs,
preconditioner);
constraints.distribute(completely_distributed_solution);
locally_relevant_solution = completely_distributed_solution;
if (time_step_n % parameters.save_interval == 0)
{
- pvtu_output.write_output(dof_handler, locally_relevant_solution,
- time_step_n, current_time);
+ pvtu_output.write_output(dof_handler,
+ locally_relevant_solution,
+ time_step_n,
+ current_time);
}
refine_mesh();
}
-template<int dim>
-void CDRProblem<dim>::refine_mesh()
+template <int dim>
+void
+CDRProblem<dim>::refine_mesh()
{
- using FunctionMap =
- std::map<types::boundary_id, const Function<dim> *>;
+ using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
- KellyErrorEstimator<dim>::estimate
- (dof_handler, QGauss<dim - 1>(fe.degree + 1), FunctionMap(),
- locally_relevant_solution, estimated_error_per_cell);
+ KellyErrorEstimator<dim>::estimate(dof_handler,
+ QGauss<dim - 1>(fe.degree + 1),
+ FunctionMap(),
+ locally_relevant_solution,
+ estimated_error_per_cell);
// This solver uses a crude refinement strategy where cells with relatively
// high errors are refined and cells with relatively low errors are
{
cell->set_refine_flag();
}
- else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <= 1e-5)
+ else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <=
+ 1e-5)
{
cell->set_coarsen_flag();
}
if (triangulation.n_levels() > parameters.max_refinement_level)
{
- for (const auto &cell :
- triangulation.cell_iterators_on_level(parameters.max_refinement_level))
+ for (const auto &cell : triangulation.cell_iterators_on_level(
+ parameters.max_refinement_level))
{
cell->clear_refine_flag();
}
// Transferring the solution between different grids is ultimately just a
// few function calls but they must be made in exactly the right order.
parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
- solution_transfer(dof_handler);
+ solution_transfer(dof_handler);
triangulation.prepare_coarsening_and_refinement();
- solution_transfer.prepare_for_coarsening_and_refinement
- (locally_relevant_solution);
+ solution_transfer.prepare_for_coarsening_and_refinement(
+ locally_relevant_solution);
triangulation.execute_coarsening_and_refinement();
setup_dofs();
// <code>locally_relevant_solution</code>, so when
// parallel::distributed::SolutionTransfer::interpolate is called it uses
// those values to populate <code>temporary</code>.
- TrilinosWrappers::MPI::Vector temporary
- (locally_owned_dofs, mpi_communicator);
+ TrilinosWrappers::MPI::Vector temporary(locally_owned_dofs, mpi_communicator);
solution_transfer.interpolate(temporary);
// After <code>temporary</code> has the correct value, this call correctly
// populates <code>completely_distributed_solution</code>, which had its
}
-template<int dim>
-void CDRProblem<dim>::run()
+template <int dim>
+void
+CDRProblem<dim>::run()
{
setup_geometry();
setup_dofs();
}
-constexpr int dim {2};
+constexpr int dim{2};
-int main(int argc, char *argv[])
+int
+main(int argc, char *argv[])
{
// One of the new features in C++11 is the <code>chrono</code> component of
// the standard library. This gives us an easy way to time the output.
auto t0 = std::chrono::high_resolution_clock::now();
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
- CDR::Parameters parameters;
+ CDR::Parameters parameters;
parameters.read_parameter_file("parameters.prm");
CDRProblem<dim> cdr_problem(parameters);
cdr_problem.run();
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
{
std::cout << "time elapsed: "
- << std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count()
- << " milliseconds."
- << std::endl;
+ << std::chrono::duration_cast<std::chrono::milliseconds>(t1 -
+ t0)
+ .count()
+ << " milliseconds." << std::endl;
}
return 0;
projects / code-gallery.git / commitdiff