If we implemented this, we would get code like this:
@code
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- local_matrix(i,j) += (k_inverse_values[q][0][0] *
- fe_values.shape_value_component(i,q,0) *
- fe_values.shape_value_component(j,q,0)
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ local_matrix(i, j) += (k_inverse_values[q][0][0] *
+ fe_values.shape_value_component(i, q, 0) *
+ fe_values.shape_value_component(j, q, 0)
+
k_inverse_values[q][0][1] *
- fe_values.shape_value_component(i,q,0) *
- fe_values.shape_value_component(j,q,1)
+ fe_values.shape_value_component(i, q, 0) *
+ fe_values.shape_value_component(j, q, 1)
+
k_inverse_values[q][1][0] *
- fe_values.shape_value_component(i,q,1) *
- fe_values.shape_value_component(j,q,0)
+ fe_values.shape_value_component(i, q, 1) *
+ fe_values.shape_value_component(j, q, 0)
+
k_inverse_values[q][1][1] *
- fe_values.shape_value_component(i,q,1) *
- fe_values.shape_value_component(j,q,1)
- ) *
- fe_values.JxW(q);
+ fe_values.shape_value_component(i, q, 1) *
+ fe_values.shape_value_component(j, q, 1)
+ ) *
+ fe_values.JxW(q);
@endcode
This is, at best, tedious, error prone, and not dimension independent. There
program we do that in the following way:
@code
- const FEValuesExtractors::Vector velocities (0);
- const FEValuesExtractors::Scalar pressure (dim);
+ const FEValuesExtractors::Vector velocities(0);
+ const FEValuesExtractors::Scalar pressure(dim);
...
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- local_matrix(i,j) += (fe_values[velocities].value (i, q) *
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ local_matrix(i,j) += (fe_values[velocities].value(i, q) *
k_inverse_values[q] *
- fe_values[velocities].value (j, q)
+ fe_values[velocities].value(j, q)
-
- fe_values[velocities].divergence (i, q) *
- fe_values[pressure].value (j, q)
+ fe_values[velocities].divergence(i, q) *
+ fe_values[pressure].value(j, q)
-
- fe_values[pressure].value (i, q) *
- fe_values[velocities].divergence (j, q)) *
+ fe_values[pressure].value(i, q) *
+ fe_values[velocities].divergence(j, q)) *
fe_values.JxW(q);
@endcode
@code
for (const auto &cell : dof_handler.active_cell_iterators())
{
- fe_values.reinit (cell);
+ fe_values.reinit(cell);
local_matrix = 0;
local_rhs = 0;
- right_hand_side.value_list (fe_values.get_quadrature_points(),
- rhs_values);
- k_inverse.value_list (fe_values.get_quadrature_points(),
- k_inverse_values);
+ right_hand_side.value_list(fe_values.get_quadrature_points(),
+ rhs_values);
+ k_inverse.value_list(fe_values.get_quadrature_points(),
+ k_inverse_values);
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
- const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q);
- const double div_phi_i_u = fe_values[velocities].divergence (i, q);
- const double phi_i_p = fe_values[pressure].value (i, q);
+ const Tensor<1,dim> phi_i_u = fe_values[velocities].value(i, q);
+ const double div_phi_i_u = fe_values[velocities].divergence(i, q);
+ const double phi_i_p = fe_values[pressure].value(i, q);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
- const Tensor<1,dim> phi_j_u = fe_values[velocities].value (j, q);
- const double div_phi_j_u = fe_values[velocities].divergence (j, q);
- const double phi_j_p = fe_values[pressure].value (j, q);
-
- local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * phi_j_u
- - div_phi_i_u * phi_j_p
- - phi_i_p * div_phi_j_u) *
- fe_values.JxW(q);
+ const Tensor<1,dim> phi_j_u = fe_values[velocities].value(j, q);
+ const double div_phi_j_u = fe_values[velocities].divergence(j, q);
+ const double phi_j_p = fe_values[pressure].value(j, q);
+
+ local_matrix(i, j) += (phi_i_u * k_inverse_values[q] * phi_j_u
+ - div_phi_i_u * phi_j_p
+ - phi_i_p * div_phi_j_u) *
+ fe_values.JxW(q);
}
local_rhs(i) += -phi_i_p *
// ...
- void SchurComplement::vmult (Vector<double> &dst,
- const Vector<double> &src) const
+ void SchurComplement::vmult(Vector<double> &dst,
+ const Vector<double> &src) const
{
- B.vmult (tmp1, src);
+ B.vmult(tmp1, src);
solver_M(M, tmp2, tmp1, preconditioner_M);
- B.Tvmult (dst, tmp2);
+ B.Tvmult(dst, tmp2);
}
};
@endcode
hand side of the first equation reads $B^TM^{-1}F-G$. This could be
implemented as follows:
@code
- Vector<double> schur_rhs (P.size());
- Vector<double> tmp (U.size());
- op_M_inv.vmult (tmp, F);
- transpose_operator(op_B).vmult (schur_rhs, tmp);
+ Vector<double> schur_rhs(P.size());
+ Vector<double> tmp(U.size());
+ op_M_inv.vmult(tmp, F);
+ transpose_operator(op_B).vmult(schur_rhs, tmp);
schur_rhs -= G;
@endcode
Again, this is a perfectly valid approach, but the fact that deal.II
@code
template <int dim>
void
-KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
- std::vector<Tensor<2,dim> > &values) const
+KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
+ std::vector<Tensor<2,dim>> &values) const
{
- AssertDimension (points.size(), values.size());
+ AssertDimension(points.size(), values.size());
- for (unsigned int p=0; p<points.size(); ++p)
+ for (unsigned int p = 0; p < points.size(); ++p)
{
- values[p].clear ();
+ values[p].clear();
const double distance_to_flowline
= std::fabs(points[p][1]-0.2*std::sin(10*points[p][0]));
- const double permeability = std::max(std::exp(-(distance_to_flowline*
+ const double permeability = std::max(std::exp(-(distance_to_flowline *
distance_to_flowline)
/ (0.1 * 0.1)),
0.001);
- for (unsigned int d=0; d<dim; ++d)
- values[p][d][d] = 1./permeability;
+ for (unsigned int d = 0; d < dim; ++d)
+ values[p][d][d] = 1./permeability;
}
}
@endcode
would look like this:
@code
template <int dim>
-class KInverse : public TensorFunction<2,dim>
+class KInverse : public TensorFunction<2, dim>
{
public:
- KInverse ();
+ KInverse();
- virtual void value_list (const std::vector<Point<dim> > &points,
- std::vector<Tensor<2,dim> > &values) const;
+ virtual void
+ value_list(const std::vector<Point<dim>> &points,
+ std::vector<Tensor<2, dim>> &values) const;
private:
- std::vector<Point<dim> > centers;
+ std::vector<Point<dim>> centers;
};
template <int dim>
-KInverse<dim>::KInverse ()
+KInverse<dim>::KInverse()
{
const unsigned int N = 40;
- centers.resize (N);
- for (unsigned int i=0; i<N; ++i)
- for (unsigned int d=0; d<dim; ++d)
- centers[i][d] = 2.*rand()/RAND_MAX-1;
+ centers.resize(N);
+ for (unsigned int i = 0; i < N; ++i)
+ for (unsigned int d = 0; d < dim; ++d)
+ centers[i][d] = (2.0 * std::rand()) / RAND_MAX - 1.0;
}
template <int dim>
void
-KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
- std::vector<Tensor<2,dim> > &values) const
+KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
+ std::vector<Tensor<2, dim>> &values) const
{
- AssertDimension (points.size(), values.size());
+ AssertDimension(points.size(), values.size());
- for (unsigned int p=0; p<points.size(); ++p)
+ for (unsigned int p = 0; p < points.size(); ++p)
{
- values[p].clear ();
+ values[p].clear();
double permeability = 0;
- for (unsigned int i=0; i<centers.size(); ++i)
+ for (unsigned int i = 0; i < centers.size(); ++i)
permeability += std::exp(-(points[p] - centers[i]).norm_square() / (0.1 * 0.1));
const double normalized_permeability
= std::max(permeability, 0.005);
- for (unsigned int d=0; d<dim; ++d)
- values[p][d][d] = 1./normalized_permeability;
+ for (unsigned int d = 0; d < dim; ++d)
+ values[p][d][d] = 1.0 / normalized_permeability;
}
}
@endcode