127 * #ifndef dealii__cdr_parameters_h
128 * #define dealii__cdr_parameters_h
130 * #include <deal.II/base/parameter_handler.h>
136 * I prefer to use the
ParameterHandler class in a slightly different way than
137 * usual: The
class Parameters creates, uses, and then destroys a
139 * of keeping it around. This is nice because now all of the
run time
140 * parameters are contained in a simple
class and it can be copied or passed
141 * around very easily.
151 *
double inner_radius;
152 *
double outer_radius;
159 *
unsigned int max_refinement_level;
181<a name=
"ann-common/include/deal.II-cdr/system_matrix.h"></a>
222 *
template <
int dim,
typename MatrixType>
230 *
MatrixType & system_matrix);
232 *
template <
int dim,
typename MatrixType>
241 *
MatrixType & system_matrix);
247<a name=
"ann-common/include/deal.II-cdr/system_matrix.templates.h"></a>
292 *
template <
int dim,
typename UpdateFunction>
302 *
auto & fe = dof_handler.get_fe();
303 *
const auto dofs_per_cell = fe.dofs_per_cell;
310 *
std::vector<types::global_dof_index> local_indices(dofs_per_cell);
312 *
for (
const auto &cell : dof_handler.active_cell_iterators())
314 *
if (cell->is_locally_owned())
318 *
cell->get_dof_indices(local_indices);
319 *
for (
unsigned int q = 0;
q < quad.size(); ++
q)
324 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
326 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
338 *
time_step / 2.0 * parameters.reaction_coefficient) *
339 *
fe_values.shape_value(i,
q) *
340 *
fe_values.shape_value(
j,
q) +
347 *
(fe_values.shape_value(i,
q) *
355 *
(fe_values.shape_grad(i,
q) *
356 *
fe_values.shape_grad(
j,
q))));
365 *
template <
int dim,
typename MatrixType>
374 *
MatrixType & system_matrix)
382 *
[&constraints, &system_matrix](
383 *
const std::vector<types::global_dof_index> &local_indices,
385 *
constraints.distribute_local_to_global(cell_matrix,
391 *
template <
int dim,
typename MatrixType>
399 *
MatrixType & system_matrix)
408 *
const std::vector<types::global_dof_index> &local_indices,
410 *
system_matrix.add(local_indices, cell_matrix);
418<a name=
"ann-common/include/deal.II-cdr/system_rhs.h"></a>
459 *
template <
int dim,
typename VectorType>
469 *
const double current_time,
476<a name=
"ann-common/include/deal.II-cdr/system_rhs.templates.h"></a>
517 *
template <
int dim,
typename VectorType>
527 *
const double current_time,
530 *
auto & fe = dof_handler.get_fe();
531 *
const auto dofs_per_cell = fe.dofs_per_cell;
533 *
(parameters.stop_time - parameters.start_time) / parameters.n_time_steps;
543 *
std::vector<types::global_dof_index> local_indices(dofs_per_cell);
545 *
const double previous_time{current_time -
time_step};
547 *
for (
const auto &cell : dof_handler.active_cell_iterators())
549 *
if (cell->is_locally_owned())
553 *
cell->get_dof_indices(local_indices);
554 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
560 *
for (
unsigned int q = 0;
q < quad.size(); ++
q)
569 *
fe_values.quadrature_point(
q));
570 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
572 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
585 *
parameters.reaction_coefficient) *
586 *
fe_values.shape_value(i,
q) *
587 *
fe_values.shape_value(
j,
q) -
594 *
(fe_values.shape_value(i,
q) *
602 *
(fe_values.shape_grad(i,
q) *
603 *
fe_values.shape_grad(
j,
q)))) *
612 *
fe_values.shape_value(i,
q));
616 *
constraints.distribute_local_to_global(
cell_rhs,
627<a name=
"ann-common/include/deal.II-cdr/write_pvtu_output.h"></a>
662 *
template <
int dim,
typename VectorType>
665 *
const VectorType & solution,
667 *
const double current_time);
678<a name=
"ann-common/include/deal.II-cdr/write_pvtu_output.templates.h"></a>
722 *
template <
int dim,
typename VectorType>
725 *
const VectorType & solution,
727 *
const double current_time)
730 *
data_out.attach_dof_handler(dof_handler);
731 *
data_out.add_data_vector(solution,
"u");
733 *
const auto &
triangulation = dof_handler.get_triangulation();
739 *
data_out.add_data_vector(
subdomain,
"subdomain");
743 *
flags.time = current_time;
751 *
data_out.set_flags(flags);
767 *
data_out.write_vtu(output);
769 *
if (this_mpi_process == 0)
772 *
for (
unsigned int i = 0;
788<a name=
"ann-common/source/parameters.cc"></a>
793 * -----------------------------------------------------------------------------
806 * -----------------------------------------------------------------------------
840 *
"Diffusion coefficient.");
844 *
"Reaction coefficient.");
849 *
"the forcing function depends on time.");
858 *
"Initial number of levels in the mesh.");
862 *
"Maximum number of levels in the mesh.");
866 *
"Finite element order.");
883 *
"Number of time steps.");
934 *
max_refinement_level =
959<a name=
"ann-common/source/system_matrix.cc"></a>
1074<a name=
"ann-common/source/system_rhs.cc"></a>
1106 *
using namespace dealii;
1117 *
const double current_time,
1129 *
const double current_time,
1141 *
const double current_time,
1153 *
const double current_time,
1159<a name=
"ann-common/source/write_pvtu_output.cc"></a>
1189 *
using namespace dealii;
1200 *
const double current_time);
1206 *
const double current_time);
1212 *
const double current_time);
1218 *
const double current_time);
1223<a name=
"ann-solver/cdr.cc"></a>
1264 *
#
include <
deal.II/distributed/grid_refinement.h>
1283 *
using namespace dealii;
1293 *
template <
int dim>
1304 *
double current_time;
1356 *
template <
int dim>
1358 *
: parameters(parameters)
1359 *
,
time_step{(parameters.stop_time - parameters.start_time) /
1360 *
parameters.n_time_steps}
1361 *
, current_time{parameters.start_time}
1365 *
, fe(parameters.fe_order)
1366 *
, quad(parameters.fe_order + 2)
1381 *
std::exp(-40 * Utilities::fixed_power<6>(p[0] - 1.5)) *
1382 *
std::exp(-40 * Utilities::fixed_power<6>(p[1]));
1385 *
,
pcout(std::cout, this_mpi_process == 0)
1391 *
template <
int dim>
1398 *
parameters.inner_radius,
1399 *
parameters.outer_radius);
1401 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1405 *
triangulation.refine_global(parameters.initial_refinement_level);
1409 *
template <
int dim>
1413 *
dof_handler.distribute_dofs(fe);
1414 *
pcout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
1416 *
locally_owned_dofs = dof_handler.locally_owned_dofs();
1419 *
constraints.
clear();
1425 *
constraints.close();
1431 *
mpi_communicator);
1435 *
template <
int dim>
1445 *
dof_handler.locally_owned_dofs(),
1449 *
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1450 *
system_matrix.reinit(locally_owned_dofs,
1452 *
mpi_communicator);
1462 *
preconditioner.initialize(system_matrix);
1466 *
template <
int dim>
1470 *
double current_time = parameters.start_time;
1494 *
solver.solve(system_matrix,
1501 *
if (
time_step_n % parameters.save_interval == 0)
1514 *
template <
int dim>
1518 *
using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
1535 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1539 *
cell->set_refine_flag();
1544 *
cell->set_coarsen_flag();
1548 *
if (
triangulation.n_levels() > parameters.max_refinement_level)
1550 *
for (
const auto &cell :
triangulation.cell_iterators_on_level(
1551 *
parameters.max_refinement_level))
1553 *
cell->clear_refine_flag();
1577 * parallel::distributed::SolutionTransfer::interpolate
is called it uses
1596 * distributed vector.
1605 *
template <
int dim>
1616 *
constexpr int dim{2};
1628 *
auto t0 = std::chrono::high_resolution_clock::now();
1632 *
parameters.read_parameter_file(
"parameters.prm");
1636 *
auto t1 = std::chrono::high_resolution_clock::now();
1639 *
std::cout <<
"time elapsed: "
1643 *
<<
" milliseconds." << std::endl;
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask={})
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
void hyper_shell(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation