Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Use extractors for local_rhs in step-22 #12988

Merged
merged 1 commit into from
Nov 24, 2021
Merged

Use extractors for local_rhs in step-22 #12988

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
3 changes: 3 additions & 0 deletions doc/news/changes/minor/20211123MarcoFeder
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
Fixed: The step-22 tutorial now uses FEValuesExtractors also for the right-hand side
<br>
(Marco Feder, 2021/11/23)
71 changes: 38 additions & 33 deletions examples/step-22/step-22.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/tensor_function.h>

#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
Expand All @@ -42,6 +43,7 @@
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>


#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
Expand Down Expand Up @@ -216,37 +218,43 @@ namespace Step22


// We implement similar functions for the right hand side which for the
// current example is simply zero:
// current example is simply zero. In the current case, what we want to
// represent is the right hand side of the velocity equation, that is, the
// body forces acting on the fluid. (The right hand side of the divergence
// equation is always zero, since we are dealing with an incompressible
// medium.) As such, the right hand side is represented by a `Tensor<1,dim>`
// object: a dim-dimensional vector, and the right tool to implement such a
// function is to derive it from `TensorFunction<1,dim>`:
template <int dim>
class RightHandSide : public Function<dim>
class RightHandSide : public TensorFunction<1, dim>
{
public:
RightHandSide()
: Function<dim>(dim + 1)
: TensorFunction<1, dim>()
{}

virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
virtual Tensor<1, dim> value(const Point<dim> &p) const override;

virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override;
virtual void value_list(const std::vector<Point<dim>> &p,
std::vector<Tensor<1, dim>> &value) const override;
};


template <int dim>
double RightHandSide<dim>::value(const Point<dim> & /*p*/,
const unsigned int /*component*/) const
Tensor<1, dim> RightHandSide<dim>::value(const Point<dim> & /*p*/) const
{
return 0;
return Tensor<1, dim>();
}


template <int dim>
void RightHandSide<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
void RightHandSide<dim>::value_list(const std::vector<Point<dim>> &vp,
std::vector<Tensor<1, dim>> &values) const
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = RightHandSide<dim>::value(p, c);
for (unsigned int c = 0; c < vp.size(); ++c)
{
values[c] = RightHandSide<dim>::value(vp[c]);
}
}


Expand Down Expand Up @@ -629,7 +637,7 @@ namespace Step22
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

const RightHandSide<dim> right_hand_side;
std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim + 1));
std::vector<Tensor<1, dim>> rhs_values(n_q_points, Tensor<1, dim>());

// Next, we need two objects that work as extractors for the FEValues
// object. Their use is explained in detail in the report on @ref
Expand Down Expand Up @@ -661,6 +669,7 @@ namespace Step22
// the outer loop.
std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
std::vector<double> div_phi_u(dofs_per_cell);
std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
std::vector<double> phi_p(dofs_per_cell);

for (const auto &cell : dof_handler.active_cell_iterators())
Expand All @@ -670,8 +679,8 @@ namespace Step22
local_preconditioner_matrix = 0;
local_rhs = 0;

right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
rhs_values);
right_hand_side.value_list(fe_values.get_quadrature_points(),
rhs_values);

for (unsigned int q = 0; q < n_q_points; ++q)
{
Expand All @@ -680,6 +689,7 @@ namespace Step22
symgrad_phi_u[k] =
fe_values[velocities].symmetric_gradient(k, q);
div_phi_u[k] = fe_values[velocities].divergence(k, q);
phi_u[k] = fe_values[velocities].value(k, q);
phi_p[k] = fe_values[pressure].value(k, q);
}

Expand Down Expand Up @@ -729,22 +739,17 @@ namespace Step22
// is overloaded for symmetric tensors, yielding the scalar
// product between the two tensors.
//
// For the right-hand side we use the fact that the shape
// functions are only non-zero in one component (because our
// elements are primitive). Instead of multiplying the tensor
// representing the dim+1 values of shape function i with the
// whole right-hand side vector, we only look at the only
// non-zero component. The function
// FiniteElement::system_to_component_index will return
// which component this shape function lives in (0=x velocity,
// 1=y velocity, 2=pressure in 2d), which we use to pick out
// the correct component of the right-hand side vector to
// multiply with.
const unsigned int component_i =
fe.system_to_component_index(i).first;
local_rhs(i) += (fe_values.shape_value(i, q) // (phi_u_i(x_q)
* rhs_values[q](component_i)) // * f(x_q))
* fe_values.JxW(q); // * dx
// For the right-hand side, we need to multiply the (vector of)
// velocity shape functions with the vector of body force
// right-hand side components, both evaluated at the current
// quadrature point. We have implemented the body forces as a
// `TensorFunction<1,dim>`, so its values at quadrature points
// are already tensors for which the application of `operator*`
// against the velocity components of the shape function results
// in the dot product, as intended.
local_rhs(i) += phi_u[i] // phi_u_i(x_q)
* rhs_values[q] // * f(x_q)
* fe_values.JxW(q); // * dx
}
}

Expand Down