<i>maximum</i> number of children a cell can have. How many children a
refined cell has was previously available as static information, but
now it depends on the actual refinement state of a cell and can be
-retrieved using the function call <code>cell-@>n_children()</code>,
+retrieved using TriaAccessor::n_children(),
a call that works equally well for both isotropic and anisotropic
refinement. A very similar situation can be found for
-faces and their subfaces: the previously available variable
-GeometryInfo<dim>::subfaces_per_face no
-longer exists; the pertinent information can now be queried using
+faces and their subfaces: the pertinent information can be queried using
GeometryInfo<dim>::max_children_per_face or <code>face->n_children()</code>,
depending on the context.
i.e. has children occupying only part of the
common face. In this case, the face
under consideration has to be a refined one, which can determine by
- asking <code>if(face->has_children())</code>. If this is true, we need to
+ asking <code>if (face->has_children())</code>. If this is true, we need to
loop over
all subfaces and get the neighbors' child behind this subface, so that we can
reinit an FEFaceValues object with the neighbor and an FESubfaceValues object
we refine the mesh and approximate the true solution better and better.
Thus, a large jump
across a given face indicates that the cell should be refined (at least)
-orthogonal to that face, whereas a small jump does not lead to this
+orthogonally to that face, whereas a small jump does not lead to this
conclusion. It is possible, of course, that the exact solution is not smooth and
that it also features a jump. In that case, however, a large jump over one face
indicates, that this face is more or less parallel to the jump and in the
// The deal.II include files have already been covered in previous examples
// and will thus not be further commented on.
-#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
-#include <deal.II/lac/vector.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/solver_richardson.h>
#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/fe/fe_values.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/numerics/data_out.h>
+#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe_dgq.h>
-#include <deal.II/lac/solver_richardson.h>
-#include <deal.II/lac/precondition_block.h>
+#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/derivative_approximation.h>
-#include <deal.II/base/timer.h>
// And this again is C++:
+#include <array>
#include <iostream>
#include <fstream>
public:
virtual void value_list(const std::vector<Point<dim>> &points,
std::vector<double> & values,
- const unsigned int component = 0) const override;
+ const unsigned int /*component*/ = 0) const override
+ {
+ (void)points;
+ Assert(values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ std::fill(values.begin(), values.end(), 0.);
+ }
};
public:
virtual void value_list(const std::vector<Point<dim>> &points,
std::vector<double> & values,
- const unsigned int component = 0) const override;
+ const unsigned int /*component*/ = 0) const override
+ {
+ Assert(values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i = 0; i < values.size(); ++i)
+ {
+ if (points[i](0) < 0.5)
+ values[i] = 1.;
+ else
+ values[i] = 0.;
+ }
+ }
};
class Beta
{
public:
+ // The flow field is chosen to be a quarter circle with counterclockwise
+ // flow direction and with the origin as midpoint for the right half of the
+ // domain with positive $x$ values, whereas the flow simply goes to the left
+ // in the left part of the domain at a velocity that matches the one coming
+ // in from the right. In the circular part the magnitude of the flow
+ // velocity is proportional to the distance from the origin. This is a
+ // difference to step-12, where the magnitude was 1 everywhere. the new
+ // definition leads to a linear variation of $\beta$ along each given face
+ // of a cell. On the other hand, the solution $u(x,y)$ is exactly the same
+ // as before.
void value_list(const std::vector<Point<dim>> &points,
- std::vector<Point<dim>> & values) const;
- };
-
-
- template <int dim>
- void RHS<dim>::value_list(const std::vector<Point<dim>> &points,
- std::vector<double> & values,
- const unsigned int) const
- {
- (void)points;
- Assert(values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- std::fill(values.begin(), values.end(), 0.);
- }
-
-
- // The flow field is chosen to be a quarter circle with counterclockwise
- // flow direction and with the origin as midpoint for the right half of the
- // domain with positive $x$ values, whereas the flow simply goes to the left
- // in the left part of the domain at a velocity that matches the one coming
- // in from the right. In the circular part the magnitude of the flow
- // velocity is proportional to the distance from the origin. This is a
- // difference to step-12, where the magnitude was 1 everywhere. the new
- // definition leads to a linear variation of $\beta$ along each given face
- // of a cell. On the other hand, the solution $u(x,y)$ is exactly the same
- // as before.
- template <int dim>
- void Beta<dim>::value_list(const std::vector<Point<dim>> &points,
- std::vector<Point<dim>> & values) const
- {
- Assert(values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i = 0; i < points.size(); ++i)
- {
- if (points[i](0) > 0)
- {
- values[i](0) = -points[i](1);
- values[i](1) = points[i](0);
- }
- else
- {
- values[i] = Point<dim>();
- values[i](0) = -points[i](1);
- }
- }
- }
-
+ std::vector<Point<dim>> & values) const
+ {
+ Assert(values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
- template <int dim>
- void BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
- std::vector<double> & values,
- const unsigned int) const
- {
- Assert(values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
+ for (unsigned int i = 0; i < points.size(); ++i)
+ {
+ if (points[i](0) > 0)
+ {
+ values[i](0) = -points[i](1);
+ values[i](1) = points[i](0);
+ }
+ else
+ {
+ values[i] = Point<dim>();
+ values[i](0) = -points[i](1);
+ }
+ }
+ }
+ };
- for (unsigned int i = 0; i < values.size(); ++i)
- {
- if (points[i](0) < 0.5)
- values[i] = 1.;
- else
- values[i] = 0.;
- }
- }
// @sect3{Class: DGTransportEquation}
};
+
// Likewise, the constructor of the class as well as the functions
// assembling the terms corresponding to cell interiors and boundary faces
// are unchanged from before. The function that assembles face terms between
{}
+
template <int dim>
void DGTransportEquation<dim>::assemble_cell_term(
const FEValues<dim> &fe_v,
{
public:
DGMethod(const bool anisotropic);
- ~DGMethod();
void run();
{}
- template <int dim>
- DGMethod<dim>::~DGMethod()
- {
- dof_handler.clear();
- }
-
template <int dim>
void DGMethod<dim>::setup_system()
// accumulate these values into vectors with
// <code>dim</code> components.
jump[face_no / 2] +=
- std::fabs(u[x] - u_neighbor[x]) * JxW[x];
+ std::abs(u[x] - u_neighbor[x]) * JxW[x];
// We also sum up the scaled weights to obtain
// the measure of the face.
area[face_no / 2] += JxW[x];
++x)
{
jump[face_no / 2] +=
- std::fabs(u[x] - u_neighbor[x]) * JxW[x];
+ std::abs(u[x] - u_neighbor[x]) * JxW[x];
area[face_no / 2] += JxW[x];
}
}
++x)
{
jump[face_no / 2] +=
- std::fabs(u[x] - u_neighbor[x]) * JxW[x];
+ std::abs(u[x] - u_neighbor[x]) * JxW[x];
area[face_no / 2] += JxW[x];
}
}
}
// Now we analyze the size of the mean jumps, which we get dividing
// the jumps by the measure of the respective faces.
- double average_jumps[dim];
- double sum_of_average_jumps = 0.;
+ std::array<double, dim> average_jumps;
+ double sum_of_average_jumps = 0.;
for (unsigned int i = 0; i < dim; ++i)
{
average_jumps[i] = jump(i) / area(i);
}
+
template <int dim>
void DGMethod<dim>::run()
{
projects / dealii.git / commitdiff