@code
hp::DoFHandler<dim> dof_handler (triangulation);
- for (typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell)
+ for (auto &cell: dof_handler.active_cell_iterators())
cell->set_active_fe_index (...);
dof_handler.distribute_dofs (fe_collection);
@endcode
conceptually very similar to how we compute hanging node constraints,
and in fact the code looks exactly the same:
@code
- ConstraintMatrix constraints;
- DoFTools::make_hanging_node_constraints (dof_handler,
- constraints);
+ AffineConstraints<double> constraints;
+ DoFTools::make_hanging_node_constraints (dof_handler, constraints);
@endcode
In other words, the DoFTools::make_hanging_node_constraints deals not
only with hanging node constraints, but also with $hp$ constraints at
update_values | update_gradients |
update_q_points | update_JxW_values);
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell!=endc; ++cell)
+ for (const auto &cell: dof_handler.active_cell_iterators())
{
hp_fe_values.reinit (cell,
cell->active_fe_index(),
Most programs built on deal.II use the DoFTools::make_sparsity_pattern
function to allocate the sparsity pattern of a matrix, and later add a few
more entries necessary to handle constrained degrees of freedom using
-ConstraintMatrix::condense. The sparsity pattern is then compressed using
+AffineConstraints::condense. The sparsity pattern is then compressed using
SparsityPattern::compress. This method is explained in step-6 and used in
most tutorial programs. In order to work, it needs an initial upper estimate
for the maximal number of nonzero entries per row, something that can be had
It turns out that the strategy presented first in step-6 to eliminate the
constraints while computing the element matrices and vectors with
-ConstraintMatrix::distribute_local_to_global is the most efficient approach
+AffineConstraints::distribute_local_to_global is the most efficient approach
also for this case. The alternative strategy to first build the matrix without
constraints and then "condensing" away constrained degrees of freedom is
considerably more expensive. It turns out that building the sparsity pattern
inhomogeneous) constraints and eliminate the matrix rows and columns to
those as well. All we have to do for this is to call the function that
interpolates the Dirichlet boundary conditions already in the setup phase in
-order to tell the ConstraintMatrix object about them, and then do the
+order to tell the AffineConstraints object about them, and then do the
transfer from local to global data on matrix and vector simultaneously. This
is exactly what we've shown in step-6.
std::vector<double> ln_k;
Table<dim, std::complex<double>> fourier_coefficients;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
std::vector<types::global_dof_index> local_dof_indices;
- typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler
- .begin_active(),
- endc = dof_handler.end();
- for (; cell != endc; ++cell)
+ for (const auto &cell : dof_handler.active_cell_iterators())
{
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
cell_matrix(i, j) +=
- (fe_values.shape_grad(i, q_point) *
- fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
+ (fe_values.shape_grad(i, q_point) * // grad phi_i(x_q)
+ fe_values.shape_grad(j, q_point) * // grad phi_j(x_q)
+ fe_values.JxW(q_point)); // dx
- cell_rhs(i) += (fe_values.shape_value(i, q_point) *
- rhs_values[q_point] * fe_values.JxW(q_point));
+ cell_rhs(i) += (fe_values.shape_value(i, q_point) * // phi_i(x_q)
+ rhs_values[q_point] * // f(x_q)
+ fe_values.JxW(q_point)); // dx
}
local_dof_indices.resize(dofs_per_cell);
// all integers, so that it what we use:
{
Vector<float> fe_degrees(triangulation.n_active_cells());
- {
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell != endc; ++cell)
- fe_degrees(cell->active_cell_index()) =
- fe_collection[cell->active_fe_index()].degree;
- }
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ fe_degrees(cell->active_cell_index()) =
+ fe_collection[cell->active_fe_index()].degree;
// With now all data vectors available -- solution, estimated errors and
// smoothness indicators, and finite element degrees --, we create a
smoothness_indicators.end()),
min_smoothness = *std::max_element(smoothness_indicators.begin(),
smoothness_indicators.end());
- {
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell != endc; ++cell)
- if (cell->refine_flag_set())
- {
- max_smoothness =
- std::max(max_smoothness,
- smoothness_indicators(cell->active_cell_index()));
- min_smoothness =
- std::min(min_smoothness,
- smoothness_indicators(cell->active_cell_index()));
- }
- }
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ if (cell->refine_flag_set())
+ {
+ max_smoothness =
+ std::max(max_smoothness,
+ smoothness_indicators(cell->active_cell_index()));
+ min_smoothness =
+ std::min(min_smoothness,
+ smoothness_indicators(cell->active_cell_index()));
+ }
const float threshold_smoothness = (max_smoothness + min_smoothness) / 2;
// With this, we can go back, loop over all cells again, and for those
// degree and in return remove the flag indicating that the cell should
// undergo bisection. For all other cells, the refinement flags remain
// untouched:
- {
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell != endc; ++cell)
- if (cell->refine_flag_set() &&
- (smoothness_indicators(cell->active_cell_index()) >
- threshold_smoothness) &&
- (cell->active_fe_index() + 1 < fe_collection.size()))
- {
- cell->clear_refine_flag();
- cell->set_active_fe_index(cell->active_fe_index() + 1);
- }
- }
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ if (cell->refine_flag_set() &&
+ (smoothness_indicators(cell->active_cell_index()) >
+ threshold_smoothness) &&
+ (cell->active_fe_index() + 1 < fe_collection.size()))
+ {
+ cell->clear_refine_flag();
+ cell->set_active_fe_index(cell->active_fe_index() + 1);
+ }
// At the end of this procedure, we then refine the mesh. During this
// process, children of cells undergoing bisection inherit their mother
{
const unsigned int dim = 2;
- static const Point<2> vertices_1[] = {
- Point<2>(-1., -1.), Point<2>(-1. / 2, -1.),
- Point<2>(0., -1.), Point<2>(+1. / 2, -1.),
- Point<2>(+1, -1.),
-
- Point<2>(-1., -1. / 2.), Point<2>(-1. / 2, -1. / 2.),
- Point<2>(0., -1. / 2.), Point<2>(+1. / 2, -1. / 2.),
- Point<2>(+1, -1. / 2.),
-
- Point<2>(-1., 0.), Point<2>(-1. / 2, 0.),
- Point<2>(+1. / 2, 0.), Point<2>(+1, 0.),
-
- Point<2>(-1., 1. / 2.), Point<2>(-1. / 2, 1. / 2.),
- Point<2>(0., 1. / 2.), Point<2>(+1. / 2, 1. / 2.),
- Point<2>(+1, 1. / 2.),
-
- Point<2>(-1., 1.), Point<2>(-1. / 2, 1.),
- Point<2>(0., 1.), Point<2>(+1. / 2, 1.),
- Point<2>(+1, 1.)};
- const unsigned int n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
- const std::vector<Point<dim>> vertices(&vertices_1[0],
- &vertices_1[n_vertices]);
- static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell] = {
- {0, 1, 5, 6},
- {1, 2, 6, 7},
- {2, 3, 7, 8},
- {3, 4, 8, 9},
- {5, 6, 10, 11},
- {8, 9, 12, 13},
- {10, 11, 14, 15},
- {12, 13, 17, 18},
- {14, 15, 19, 20},
- {15, 16, 20, 21},
- {16, 17, 21, 22},
- {17, 18, 22, 23}};
- const unsigned int n_cells =
- sizeof(cell_vertices) / sizeof(cell_vertices[0]);
+ const std::vector<Point<2>> vertices = {
+ {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0}, //
+ {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5}, //
+ {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0}, //
+ {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5}, //
+ {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}};
+
+ const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
+ cell_vertices = {{0, 1, 5, 6},
+ {1, 2, 6, 7},
+ {2, 3, 7, 8},
+ {3, 4, 8, 9},
+ {5, 6, 10, 11},
+ {8, 9, 12, 13},
+ {10, 11, 14, 15},
+ {12, 13, 17, 18},
+ {14, 15, 19, 20},
+ {15, 16, 20, 21},
+ {16, 17, 21, 22},
+ {17, 18, 22, 23}};
+
+ const unsigned int n_cells = cell_vertices.size();
std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
for (unsigned int i = 0; i < n_cells; ++i)
Vector<double> local_dof_values;
// Then here is the loop:
- typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler
- .begin_active(),
- endc = dof_handler.end();
- for (; cell != endc; ++cell)
+ for (const auto &cell : dof_handler.active_cell_iterators())
{
// Inside the loop, we first need to get the values of the local
// degrees of freedom (which we put into the
// We have to calculate the logarithms of absolute
// values of coefficients and use it in linear regression fit to
// obtain $\mu$.
- for (unsigned int f = 0; f < res.second.size(); f++)
- res.second[f] = std::log(res.second[f]);
+ for (double &f : res.second)
+ f = std::log(f);
std::pair<double, double> fit =
FESeries::linear_regression(ln_k, res.second);