$u_h = \sum_{i=1}^N U_i \varphi_i$. So using this we can give an expression for
the discrete Jacobian and the residual:
@f{align*}{
- A_{i,j} = \bigl( F'(u_h^n) \bigr)_{i,j}
+ A_{i,j} = \bigl( F'(u_h^n) \bigr)_{i,j}
&=
\int_\Omega \nabla\varphi_i \cdot \nabla \varphi_j \,\mathrm{d} x
-
\int_\Omega \varphi_i \, \exp( u_h ) \varphi_j \,\mathrm{d} x,\\
- b_{i} = \bigl( F(u_h^n) \bigr)_{i}
+ b_{i} = \bigl( F(u_h^n) \bigr)_{i}
&=
\int_\Omega \nabla\varphi_i \cdot \nabla u_h^n \,\mathrm{d} x
-
{
cell_matrix = 0.0;
cell_rhs = 0.0;
-
+
fe_values.reinit(cell);
-
+
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 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->get_dof_indices(local_dof_indices);
-
+
constraints.distribute_local_to_global(cell_matrix, cell_rhs,
local_dof_indices,
system_matrix, system_rhs);
-
+
}
-
+
}
@endcode
GelfandProblem<dim>::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]);
}
-
+
// Define the actual preconditioner
...
-
+
// Solve the linear system
...
}
{
const unsigned int n_cells = this->data->n_cell_batches();
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 q = 0; q < phi.n_q_points; ++q)
nonlinear_values(cell, q) = std::exp(phi.get_value(q));
}
JacobianOperator();
- virtual void
- clear() override;
+ virtual void clear() override;
- void
- evaluate_newton_step(
+ void evaluate_newton_step(
const LinearAlgebra::distributed::Vector<number> &newton_step);
- virtual void
- compute_diagonal() override;
+ virtual void compute_diagonal() override;
private:
- virtual void
- apply_add(
+ virtual void apply_add(
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src) const override;
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;
- void
- local_compute_diagonal(FECellIntegrator &integrator) const;
+ void local_compute_diagonal(FECellIntegrator &integrator) const;
Table<2, VectorizedArray<number>> nonlinear_values;
};
// the nonlinearity and call the <code>clear()</code> function of the base
// class.
template <int dim, int fe_degree, typename number>
- void
- JacobianOperator<dim, fe_degree, number>::clear()
+ void JacobianOperator<dim, fe_degree, number>::clear()
{
nonlinear_values.reinit(0, 0);
MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>::
// This skips all evaluations of the nonlinearity in each call of the
// <code>vmult()</code> function.
template <int dim, int fe_degree, typename number>
- void
- JacobianOperator<dim, fe_degree, number>::evaluate_newton_step(
+ 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();
// to perform the cell integration and distribute the local contributions into
// the global vector <code> dst</code>.
template <int dim, int fe_degree, typename number>
- void
- JacobianOperator<dim, fe_degree, number>::local_apply(
+ void JacobianOperator<dim, fe_degree, number>::local_apply(
const MatrixFree<dim, number> & data,
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src,
// Next we use MatrixFree::cell_loop() to perform the actual loop over all
// cells computing the cell contribution to the matrix-vector product.
template <int dim, int fe_degree, typename number>
- void
- JacobianOperator<dim, fe_degree, number>::apply_add(
+ void JacobianOperator<dim, fe_degree, number>::apply_add(
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src) const
{
// values from a input vector or distribute any local results to an output
// vector. Instead the only input argument is the used FEEvaluation object.
template <int dim, int fe_degree, typename number>
- void
- JacobianOperator<dim, fe_degree, number>::local_compute_diagonal(
+ void JacobianOperator<dim, fe_degree, number>::local_compute_diagonal(
FECellIntegrator &phi) const
{
AssertDimension(nonlinear_values.size(0),
// diagonal and invert the elements by hand. Note, that during this loop we
// catch the constrained DOFs and set them manually to one.
template <int dim, int fe_degree, typename number>
- void
- JacobianOperator<dim, fe_degree, number>::compute_diagonal()
+ void JacobianOperator<dim, fe_degree, number>::compute_diagonal()
{
this->inverse_diagonal_entries.reset(
new DiagonalMatrix<LinearAlgebra::distributed::Vector<number>>());
public:
GelfandProblem();
- void
- run();
+ void run();
private:
- void
- make_grid();
+ void make_grid();
- void
- setup_system();
+ void setup_system();
- void
- evaluate_residual(
+ void evaluate_residual(
LinearAlgebra::distributed::Vector<double> & dst,
const LinearAlgebra::distributed::Vector<double> &src) const;
- void
- local_evaluate_residual(
+ void local_evaluate_residual(
const MatrixFree<dim, double> & data,
LinearAlgebra::distributed::Vector<double> & dst,
const LinearAlgebra::distributed::Vector<double> &src,
const std::pair<unsigned int, unsigned int> & cell_range) const;
- void
- assemble_rhs();
+ void assemble_rhs();
- double
- compute_residual(const double alpha);
+ double compute_residual(const double alpha);
- void
- compute_update();
+ void compute_update();
- void
- solve();
+ void solve();
- double
- compute_solution_norm() const;
+ double compute_solution_norm() const;
- void
- output_results(const unsigned int cycle) const;
+ void output_results(const unsigned int cycle) const;
// For the parallel computation we define a
// class and also assign a SphericalManifold for the boundary. Finally, we
// refine the initial mesh 3 - <code>dim</code> times globally.
template <int dim, int fe_degree>
- void
- GelfandProblem<dim, fe_degree>::make_grid()
+ void GelfandProblem<dim, fe_degree>::make_grid()
{
TimerOutput::Scope t(computing_timer, "make grid");
// Note how we can use the same MatrixFree object twice, for the
// <code>JacobianOperator</code> and the multigrid preconditioner.
template <int dim, int fe_degree>
- void
- GelfandProblem<dim, fe_degree>::setup_system()
+ void GelfandProblem<dim, fe_degree>::setup_system()
{
TimerOutput::Scope t(computing_timer, "setup system");
// related data exchange, since all the bookkeeping is done by the
// MatrixFree::cell_loop().
template <int dim, int fe_degree>
- void
- GelfandProblem<dim, fe_degree>::evaluate_residual(
+ void GelfandProblem<dim, fe_degree>::evaluate_residual(
LinearAlgebra::distributed::Vector<double> & dst,
const LinearAlgebra::distributed::Vector<double> &src) const
{
// FEEvaluation::read_dof_values_plain() and FEEvaluation::evaluate(), since
// the input vector might have constrained DOFs.
template <int dim, int fe_degree>
- void
- GelfandProblem<dim, fe_degree>::local_evaluate_residual(
+ void GelfandProblem<dim, fe_degree>::local_evaluate_residual(
const MatrixFree<dim, double> & data,
LinearAlgebra::distributed::Vector<double> & dst,
const LinearAlgebra::distributed::Vector<double> &src,
// Experiences show that using the FEEvaluation class is much faster than a
// classical implementation with FEValues and co.
template <int dim, int fe_degree>
- void
- GelfandProblem<dim, fe_degree>::assemble_rhs()
+ void GelfandProblem<dim, fe_degree>::assemble_rhs()
{
TimerOutput::Scope t(computing_timer, "assemble right hand side");
// use a damped version $\alpha<1$ until the Newton step is good enough and
// the full Newton step can be performed. This was also discussed in step-15.
template <int dim, int fe_degree>
- double
- GelfandProblem<dim, fe_degree>::compute_residual(const double alpha)
+ double GelfandProblem<dim, fe_degree>::compute_residual(const double alpha)
{
TimerOutput::Scope t(computing_timer, "compute residual");
// preconditioner. For this we first set up the PreconditionMG object with a
// Chebyshev smoother like we did in step-37.
template <int dim, int fe_degree>
- void
- GelfandProblem<dim, fe_degree>::compute_update()
+ void GelfandProblem<dim, fe_degree>::compute_update()
{
TimerOutput::Scope t(computing_timer, "compute update");
// Now we implement the actual Newton solver for the nonlinear problem.
template <int dim, int fe_degree>
- void
- GelfandProblem<dim, fe_degree>::solve()
+ void GelfandProblem<dim, fe_degree>::solve()
{
TimerOutput::Scope t(computing_timer, "solve");
// VectorTools::integrate_difference(). In the end we gather all computations
// from all MPI ranks and return the norm.
template <int dim, int fe_degree>
- double
- GelfandProblem<dim, fe_degree>::compute_solution_norm() const
+ double GelfandProblem<dim, fe_degree>::compute_solution_norm() const
{
solution.update_ghost_values();
// about the system specifications and the finite element space we use. The
// problem is solved several times on a successively refined mesh.
template <int dim, int fe_degree>
- void
- GelfandProblem<dim, fe_degree>::run()
+ void GelfandProblem<dim, fe_degree>::run()
{
{
const unsigned int n_ranks =
// create an object of the <code>GelfandProblem</code> class and call the run
// function. Exemplarily we solve the problem once in 2D and once in 3D each
// with fourth-order Lagrangian finite elements.
-int
-main(int argc, char *argv[])
+int main(int argc, char *argv[])
{
try
{