FEValuesBase::get_function_gradients(). The <code>assemble_system()</code>
function would then looks like:
@code
-template<int dim>
-void
-GelfandProblem<dim>::assemble_system()
+template <int dim>
+void GelfandProblem<dim>::assemble_system()
{
system_matrix = 0;
system_rhs = 0;
- const QGauss<dim> quadrature_formula(fe.degree+1);
- FEValues<dim> fe_values(fe, quadrature_formula,
+ const QGauss<dim> quadrature_formula(fe.degree + 1);
+ FEValues<dim> fe_values(fe,
+ quadrature_formula,
update_values | update_gradients | update_JxW_values);
const unsigned int n_q_points = fe_values.n_quadrature_points;
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
- std::vector<Tensor<1,dim> > newton_step_gradients(n_q_points);
+ std::vector<Tensor<1, dim>> newton_step_gradients(n_q_points);
std::vector<double> newton_step_values(n_q_points);
- for(const auto &cell : dof_handler.active_cell_iterators())
- {
- cell_matrix = 0.0;
- cell_rhs = 0.0;
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ cell_matrix = 0.0;
+ cell_rhs = 0.0;
- fe_values.reinit(cell);
+ fe_values.reinit(cell);
- fe_values.get_function_values(solution, newton_step_values);
- fe_values.get_function_gradients(solution, newton_step_gradients);
+ fe_values.get_function_values(solution, newton_step_values);
+ fe_values.get_function_gradients(solution, newton_step_gradients);
- for(unsigned int q=0; q<n_q_points; ++q)
- {
- const double nonlinearity = std::exp(newton_step_values[q]);
- const double dx = fe_values.JxW(q);
- for(unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double phi_i = fe_values.shape_value(i,q);
- const Tensor<1,dim> grad_phi_i = fe_values.shape_grad(i,q);
- for(unsigned int j=0; j<dofs_per_cell; ++j)
+ for (unsigned int q = 0; q < n_q_points; ++q)
{
- const double phi_j = fe_values.shape_value(j,q);
- const Tensor<1,dim> grad_phi_j = fe_values.shape_grad(j,q);
-
- cell_matrix(i,j) += ( grad_phi_i*grad_phi_j - phi_i*nonlinearity*phi_j ) * dx;
+ const double nonlinearity = std::exp(newton_step_values[q]);
+ const double dx = fe_values.JxW(q);
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ const double phi_i = fe_values.shape_value(i, q);
+ const Tensor<1, dim> grad_phi_i = fe_values.shape_grad(i, q);
+
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ {
+ const double phi_j = fe_values.shape_value(j, q);
+ const Tensor<1, dim> grad_phi_j = fe_values.shape_grad(j, q);
+
+ cell_matrix(i, j) +=
+ (grad_phi_i * grad_phi_j - phi_i * nonlinearity * phi_j) *
+ dx;
+ }
+
+ cell_rhs(i) += (-grad_phi_i * newton_step_gradients[q] +
+ phi_i * newton_step_values[q]) *
+ dx;
+ }
}
- cell_rhs(i) += ( -grad_phi_i*newton_step_gradients[q] + phi_i*newton_step_values[q] ) * dx;
+ cell->get_dof_indices(local_dof_indices);
- }
+ constraints.distribute_local_to_global(
+ cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
-
- cell->get_dof_indices(local_dof_indices);
-
- constraints.distribute_local_to_global(cell_matrix, cell_rhs,
- local_dof_indices,
- system_matrix, system_rhs);
-
- }
-
}
@endcode
of the triangulation. In the code this task will be done by the function
MGTransferMatrixFree::interpolate_to_mg():
@code
-void
-GelfandProblem<dim>::compute_update()
+template <int dim, int fe_degree>
+void GelfandProblem<dim, fe_degree>::compute_update()
{
+ TimerOutput::Scope t(computing_timer, "compute update");
+
solution.update_ghost_values();
system_matrix.evaluate_newton_step(solution);
- MGTransferMatrixFree<dim,float> mg_transfer(mg_constrained_dofs);
-
mg_transfer.interpolate_to_mg(dof_handler, mg_solution, solution);
+
// Set up options for the multilevel preconditioner
- for(unsigned int level=0; level<triangulation.n_levels()-1; ++level)
- {
- mg_matrices[level].evaluate_newton_step(mg_solution[level]);
- }
+ // ...
+
+ for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
+ {
+ mg_matrices[level].evaluate_newton_step(mg_solution[level]);
+ }
// Define the actual preconditioner
- ...
+ // ...
// Solve the linear system
- ...
+ // ...
}
@endcode
function. The idea is to use an FEEvaluation object to evaluate the Newton step
and store the expression in a table for all cells and all quadrature points:
@code
-void
-JacobianOperator<dim,fe_degree,number>::evaluate_newton_step(const LinearAlgebra::distributed::Vector<number> &src)
+template <int dim, int fe_degree, typename number>
+void JacobianOperator<dim, fe_degree, number>::evaluate_newton_step(
+ const LinearAlgebra::distributed::Vector<number> &newton_step)
{
const unsigned int n_cells = this->data->n_cell_batches();
- FEEvaluation<dim, fe_degree, fe_degree+1, 1, number> phi(*this->data);
+
+ FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*this->data);
nonlinear_values.reinit(n_cells, phi.n_q_points);
- for(unsigned int cell=0; cell<n_cells; ++cell)
- {
- phi.reinit(cell);
- phi.read_dof_values_plain(src);
- phi.evaluate(true, false);
+ for (unsigned int cell = 0; cell < n_cells; ++cell)
+ {
+ phi.reinit(cell);
+ phi.read_dof_values_plain(newton_step);
+ phi.evaluate(EvaluationFlags::values);
- for (unsigned int q = 0; q < phi.n_q_points; ++q)
- nonlinear_values(cell, q) = std::exp(phi.get_value(q));
- }
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ {
+ nonlinear_values(cell, q) = std::exp(phi.get_value(q));
+ }
+ }
}
@endcode