RefinementWeightedKelly<dim>::refine_grid ()
{
// First compute some residual based error indicators for all cells by a
- // method already implemented in the library. What exactly is computed
- // can be read in the documentation of that class.
+ // method already implemented in the library. What exactly we compute
+ // here is described in more detail in the documentation of that class.
Vector<float> estimated_error_per_cell (this->triangulation->n_active_cells());
KellyErrorEstimator<dim>::estimate (this->dof_handler,
*this->face_quadrature,
this->solution,
estimated_error_per_cell);
- // Now we are going to weight these indicators by the value of the
- // function given to the constructor:
+ // Next weigh each entry in the vector of indicators by the value of the
+ // function given to the constructor, evaluated at the cell center. We need
+ // to write the result into the vector entry that corresponds to the current
+ // cell, which we can obtain by asking the cell what its index among all
+ // active cells is using CellAccessor::active_cell_index(). (In reality,
+ // this index is zero for the first cell we handle in the loop, one for the
+ // second cell, etc., and we could as well just keep track of this index
+ // using an integer counter; but using CellAccessor::active_cell_index()
+ // makes this more explicit.)
typename DoFHandler<dim>::active_cell_iterator
cell = this->dof_handler.begin_active(),
endc = this->dof_handler.end();
- for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
- estimated_error_per_cell(cell_index)
+ for (; cell!=endc; ++cell)
+ estimated_error_per_cell(cell->active_cell_index())
*= weighting_function->value (cell->center());
GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
- for (unsigned int index=0; cell!=endc; ++cell, ++index)
+ for (; cell!=endc; ++cell)
if (cell->is_locally_owned() )
{
// On these cells, add up the stresses over all quadrature
.old_stress;
// ...then write the norm of the average to their destination:
- norm_of_stress(index)
+ norm_of_stress(cell->active_cell_index())
= (accumulated_stress /
quadrature_formula.size()).norm();
}
// elements would not appear in the output file, that we would find out
// by looking at the graphical output:
else
- norm_of_stress(index) = -1e+20;
+ norm_of_stress(cell->active_cell_index()) = -1e+20;
}
// Finally attach this vector as well to be treated for output:
data_out.add_data_vector (norm_of_stress, "norm_of_stress");
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
- for (unsigned int index=0; cell!=endc; ++cell, ++index)
- fe_degrees(index)
+ for (; cell!=endc; ++cell)
+ fe_degrees(cell->active_cell_index())
= fe_collection[cell->active_fe_index()].degree;
}
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
- for (unsigned int index=0; cell!=endc; ++cell, ++index)
+ for (; cell!=endc; ++cell)
if (cell->refine_flag_set())
{
max_smoothness = std::max (max_smoothness,
- smoothness_indicators(index));
+ smoothness_indicators(cell->active_cell_index()));
min_smoothness = std::min (min_smoothness,
- smoothness_indicators(index));
+ smoothness_indicators(cell->active_cell_index()));
}
}
const float threshold_smoothness = (max_smoothness + min_smoothness) / 2;
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
- for (unsigned int index=0; cell!=endc; ++cell, ++index)
+ for (; cell!=endc; ++cell)
if (cell->refine_flag_set()
&&
- (smoothness_indicators(index) > threshold_smoothness)
+ (smoothness_indicators(cell->active_cell_index()) > threshold_smoothness)
&&
(cell->active_fe_index()+1 < fe_collection.size()))
{
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
- for (unsigned int index=0; cell!=endc; ++cell, ++index)
+ for (; cell!=endc; ++cell)
{
// Inside the loop, we first need to get the values of the local
// degrees of freedom (which we put into the
// The final step is to compute the Sobolev index $s=\mu-\frac d2$ and
// store it in the vector of estimated values for each cell:
- smoothness_indicators(index) = mu - 1.*dim/2;
+ smoothness_indicators(cell->active_cell_index()) = mu - 1.*dim/2;
}
}
}
cell = triangulation.begin_active(),
endc = triangulation.end();
- for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
- if (error_indicators(cell_index) > refine_threshold)
+ for (; cell!=endc; ++cell)
+ if (error_indicators(cell->active_cell_index()) > refine_threshold)
cell->set_refine_flag ();
- else if (error_indicators(cell_index) < coarsen_threshold)
+ else if (error_indicators(cell->active_cell_index()) < coarsen_threshold)
cell->set_coarsen_flag ();
SolutionTransfer<dim> soltrans(dof_handler);
// more. The structure of these nested conditions is much the same as we
// encountered when assembling interface terms in
// <code>assemble_system</code>.
- {
- unsigned int cell_index = 0;
- for (typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell, ++cell_index)
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (cell_is_in_solid_domain (cell))
- {
- if ((cell->at_boundary(f) == false)
+ for (typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active();
+ cell != dof_handler.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell_is_in_solid_domain (cell))
+ {
+ if ((cell->at_boundary(f) == false)
+ &&
+ (((cell->neighbor(f)->level() == cell->level())
&&
- (((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == false)
- &&
- cell_is_in_fluid_domain (cell->neighbor(f)))
- ||
- ((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == true)
- &&
- (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
- (f, 0))))
- ||
- (cell->neighbor_is_coarser(f)
- &&
- cell_is_in_fluid_domain(cell->neighbor(f)))
- ))
- estimated_error_per_cell(cell_index) = 0;
- }
- else
- {
- if ((cell->at_boundary(f) == false)
+ (cell->neighbor(f)->has_children() == false)
&&
- (((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == false)
- &&
- cell_is_in_solid_domain (cell->neighbor(f)))
- ||
- ((cell->neighbor(f)->level() == cell->level())
- &&
- (cell->neighbor(f)->has_children() == true)
- &&
- (cell_is_in_solid_domain (cell->neighbor_child_on_subface
- (f, 0))))
- ||
- (cell->neighbor_is_coarser(f)
- &&
- cell_is_in_solid_domain(cell->neighbor(f)))
- ))
- estimated_error_per_cell(cell_index) = 0;
- }
- }
+ cell_is_in_fluid_domain (cell->neighbor(f)))
+ ||
+ ((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == true)
+ &&
+ (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
+ (f, 0))))
+ ||
+ (cell->neighbor_is_coarser(f)
+ &&
+ cell_is_in_fluid_domain(cell->neighbor(f)))
+ ))
+ estimated_error_per_cell(cell->active_cell_index()) = 0;
+ }
+ else
+ {
+ if ((cell->at_boundary(f) == false)
+ &&
+ (((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == false)
+ &&
+ cell_is_in_solid_domain (cell->neighbor(f)))
+ ||
+ ((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == true)
+ &&
+ (cell_is_in_solid_domain (cell->neighbor_child_on_subface
+ (f, 0))))
+ ||
+ (cell->neighbor_is_coarser(f)
+ &&
+ cell_is_in_solid_domain(cell->neighbor(f)))
+ ))
+ estimated_error_per_cell(cell->active_cell_index()) = 0;
+ }
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
projects / dealii.git / commitdiff