// argument denotes a mask
// selecting which components of
// vector valued finite elements we
- // want to be considered. Since we
+ // want to be considered. This sort
+ // of information is encoded using
+ // the ComponentMask class (see also
+ // @ref GlossComponentMask). Since we
// have a scalar finite element
- // anyway, this mask consists of
- // only one entry, and its value
- // must be <code>true</code>.
+ // anyway, this mask in reality should
+ // have only one entry with a
+ // <code>true</code> value. However,
+ // the ComponentMask class has
+ // semantics that allow it to
+ // represents a mask of indefinite
+ // size whose every element equals
+ // <code>true</code> when one just
+ // default constructs such an object,
+ // so this is what we'll do here.
std::vector<bool> boundary_dofs (dof_handler.n_dofs(), false);
- DoFTools::extract_boundary_dofs (dof_handler, std::vector<bool>(1,true),
+ DoFTools::extract_boundary_dofs (dof_handler,
+ ComponentMask(),
boundary_dofs);
// Now first for the generation of
// belong to the boundary and then loop
// over all of those and set the residual
// entry to zero. This happens in the
- // following lines:
+ // following lines which we have already
+ // seen used in step-11:
hanging_node_constraints.condense (residual);
std::vector<bool> boundary_dofs (dof_handler.n_dofs());
DoFTools::extract_boundary_dofs (dof_handler,
- std::vector<bool>(1,true),
+ ComponentMask(),
boundary_dofs);
for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
if (boundary_dofs[i] == true)
typename FunctionMap<dim>::type(),
localized_solution,
local_error_per_cell,
- std::vector<bool>(),
+ ComponentMask(),
0,
multithread_info.n_default_threads,
this_mpi_process);
// restrict vertical motion, it
// has only its last component
// set:
- std::vector<bool> z_component (dim, false);
- z_component[dim-1] = true;
+ FEValuesExtractors::Scalar z_component (dim-1);
std::map<unsigned int,double> boundary_values;
VectorTools::
interpolate_boundary_values (dof_handler,
IncrementalBoundaryValues<dim>(present_time,
present_timestep),
boundary_values,
- z_component);
+ fe.component_mask(z_component));
PETScWrappers::MPI::Vector tmp (mpi_communicator, dof_handler.n_dofs(),
n_local_dofs);
typename FunctionMap<dim>::type(),
incremental_displacement,
error_per_cell,
- std::vector<bool>(),
+ ComponentMask(),
0,
multithread_info.n_default_threads,
this_mpi_process);
// Now comes the implementation of
// Dirichlet boundary conditions, which
// should be evident after the discussion
- // in the introduction. All that changed is
- // that the function already appears in the
- // setup functions, whereas we were used to
- // see it in some assembly routine. Further
- // down below where we set up the mesh, we
- // will associate the top boundary where we
- // impose Dirichlet boundary conditions
- // with boundary indicator 1. We will have
- // to pass this boundary indicator as
- // second argument to the function below
- // interpolating boundary values. There is
- // one more thing, though. The function
- // describing the Dirichlet conditions was
- // defined for all components, both
- // velocity and pressure. However, the
- // Dirichlet conditions are to be set for
- // the velocity only. To this end, we use
- // a <code>component_mask</code> that
- // filters out the pressure component, so
- // that the condensation is performed on
- // velocity degrees of freedom only. Since
- // we use adaptively refined grids the
+ // in the introduction. All that changed
+ // is that the function already appears
+ // in the setup functions, whereas we
+ // were used to see it in some assembly
+ // routine. Further down below where we
+ // set up the mesh, we will associate the
+ // top boundary where we impose Dirichlet
+ // boundary conditions with boundary
+ // indicator 1. We will have to pass
+ // this boundary indicator as second
+ // argument to the function below
+ // interpolating boundary values. There
+ // is one more thing, though. The
+ // function describing the Dirichlet
+ // conditions was defined for all
+ // components, both velocity and
+ // pressure. However, the Dirichlet
+ // conditions are to be set for the
+ // velocity only. To this end, we use a
+ // ComponentMask that only selects the
+ // velocity components. The component
+ // mask is obtained from the finite
+ // element by specifying the particular
+ // components we want. Since we use
+ // adaptively refined grids the
// constraint matrix needs to be first
// filled with hanging node constraints
- // generated from the DoF handler. Note the
- // order of the two functions — we
- // first compute the hanging node
+ // generated from the DoF handler. Note
+ // the order of the two functions —
+ // we first compute the hanging node
// constraints, and then insert the
// boundary values into the constraint
- // matrix. This makes sure that we respect
- // H<sup>1</sup> conformity on boundaries
- // with hanging nodes (in three space
- // dimensions), where the hanging node
- // needs to dominate the Dirichlet boundary
- // values.
+ // matrix. This makes sure that we
+ // respect H<sup>1</sup> conformity on
+ // boundaries with hanging nodes (in
+ // three space dimensions), where the
+ // hanging node needs to dominate the
+ // Dirichlet boundary values.
{
constraints.clear ();
- std::vector<bool> component_mask (dim+1, true);
- component_mask[dim] = false;
+
+ FEValuesExtractors::Vector velocities(0);
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
VectorTools::interpolate_boundary_values (dof_handler,
1,
BoundaryValues<dim>(),
constraints,
- component_mask);
+ fe.component_mask(velocities));
}
constraints.close ();
// step-6, with the exception that we base
// the refinement only on the change in
// pressure, i.e., we call the Kelly error
- // estimator with a mask
- // object. Additionally, we do not coarsen
+ // estimator with a mask object of type
+ // ComponentMask that selects the single
+ // scalar component for the pressure that
+ // we are interested in (we get such a mask
+ // from the finite element class by
+ // specifying the component we
+ // want). Additionally, we do not coarsen
// the grid again:
template <int dim>
void
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
- std::vector<bool> component_mask (dim+1, false);
- component_mask[dim] = true;
+ FEValuesExtractors::Scalar pressure(dim);
KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1>(degree+1),
typename FunctionMap<dim>::type(),
solution,
estimated_error_per_cell,
- component_mask);
+ fe.component_mask(pressure));
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
(new TrilinosWrappers::PreconditionAMG());
std::vector<std::vector<bool> > constant_modes;
- std::vector<bool> velocity_components (dim+1,true);
- velocity_components[dim] = false;
- DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components,
+ FEValuesExtractors::Vector velocity_components(0);
+ DoFTools::extract_constant_modes (stokes_dof_handler,
+ stokes_fe.component_mask(velocity_components),
constant_modes);
TrilinosWrappers::PreconditionAMG::AdditionalData amg_data;
amg_data.constant_modes = constant_modes;
DoFTools::make_hanging_node_constraints (stokes_dof_handler,
stokes_constraints);
- std::vector<bool> velocity_mask (dim+1, true);
- velocity_mask[dim] = false;
+ FEValuesExtractors::Vector velocity_components(0);
VectorTools::interpolate_boundary_values (stokes_dof_handler,
0,
ZeroFunction<dim>(dim+1),
stokes_constraints,
- velocity_mask);
+ stokes_fe.component_mask(velocity_components));
std::set<types::boundary_id> no_normal_flux_boundaries;
no_normal_flux_boundaries.insert (1);
assemble_stokes_preconditioner ();
std::vector<std::vector<bool> > constant_modes;
- std::vector<bool> velocity_components (dim+1,true);
- velocity_components[dim] = false;
- DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components,
+ FEValuesExtractors::Vector velocity_components(0);
+ DoFTools::extract_constant_modes (stokes_dof_handler,
+ stokes_fe.component_mask(velocity_components),
constant_modes);
Mp_preconditioner.reset (new TrilinosWrappers::PreconditionJacobi());
typename FunctionMap<dim>::type(),
temperature_solution,
estimated_error_per_cell,
- std::vector<bool>(),
+ ComponentMask(),
0,
0,
triangulation.locally_owned_subdomain());
constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior);
}
+
+
template <int dim>
void PlasticityContactProblem<dim>::dirichlet_constraints ()
{
constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs);
constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes);
- std::vector<bool> component_mask (dim, true);
- component_mask[0] = true;
- component_mask[1] = true;
- component_mask[2] = true;
+ // interpolate all components of the solution
VectorTools::interpolate_boundary_values (dof_handler,
6,
EquationData::BoundaryValues<dim>(),
constraints_dirichlet_hanging_nodes,
- component_mask);
+ ComponentMask());
- component_mask[0] = true;
- component_mask[1] = true;
- component_mask[2] = false;
+ // interpolate x- and y-components of the
+ // solution (this is a bit mask, so apply
+ // operator| )
+ FEValuesExtractors::Scalar x_displacement(0);
+ FEValuesExtractors::Scalar y_displacement(1);
VectorTools::interpolate_boundary_values (dof_handler,
8,
EquationData::BoundaryValues<dim>(),
constraints_dirichlet_hanging_nodes,
- component_mask);
+ (fe.component_mask(x_displacement)
+ |
+ fe.component_mask(y_displacement)));
constraints_dirichlet_hanging_nodes.close ();
}
+
+
template <int dim>
void PlasticityContactProblem<dim>::solve ()
{
solution = distributed_solution;
}
+
+
template <int dim>
void PlasticityContactProblem<dim>::solve_newton ()
{
Timer t;
std::vector<std::vector<bool> > constant_modes;
- std::vector<bool> components (dim,true);
- components[dim] = false;
- DoFTools::extract_constant_modes (dof_handler, components,
+ DoFTools::extract_constant_modes (dof_handler,
+ ComponentMask(),
constant_modes);
additional_data.elliptic = true;
pcout<< "%%%%%% Rechenzeit output = " << run_time[4] <<std::endl;
}
+
+
template <int dim>
void PlasticityContactProblem<dim>::refine_grid ()
{
}
+
+
template <int dim>
void PlasticityContactProblem<dim>::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const
{
}
}
+
+
template <int dim>
void PlasticityContactProblem<dim>::output_results (const std::string& title) const
{
move_mesh (tmp);
}
+
+
template <int dim>
void PlasticityContactProblem<dim>::run ()
{
pcout << "Solving problem in " << dim << " space dimensions." << std::endl;
-
+
Timer t;
run_time.resize (8);
{
darcy_preconditioner_constraints.clear ();
- std::vector<bool> component_mask (dim+1, false);
- component_mask[dim] = true;
-
+ FEValuesExtractors::Scalar pressure(dim);
DoFTools::make_hanging_node_constraints (darcy_dof_handler, darcy_preconditioner_constraints);
- DoFTools::make_zero_boundary_constraints (darcy_dof_handler, darcy_preconditioner_constraints, component_mask);
+ DoFTools::make_zero_boundary_constraints (darcy_dof_handler, darcy_preconditioner_constraints,
+ darcy_fe.component_mask(pressure));
darcy_preconditioner_constraints.close ();
}
// +z face has an the applied pressure but
// is also constrained in the x- and
// y-directions.
+ //
+ // In the following, we will have to tell
+ // the function interpolation boundary
+ // values which components of the
+ // solution vector should be constrained
+ // (i.e., whether it's the x-, y-,
+ // z-displacements or combinations
+ // thereof). This is done using
+ // ComponentMask objects (see @ref
+ // GlossComponentMask) which we can get
+ // from the finite element if we provide
+ // it with an extractor object for the
+ // component we wish to select. To this
+ // end we first set up such extractor
+ // objects and later use it when
+ // generating the relevant component
+ // masks:
+ const FEValuesExtractors::Scalar x_displacement(0);
+ const FEValuesExtractors::Scalar y_displacement(1);
+ const FEValuesExtractors::Scalar z_displacement(2);
+
{
const int boundary_id = 0;
- std::vector<bool> components(n_components, false);
- components[0] = true;
-
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ fe.component_mask(x_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ fe.component_mask(x_displacement));
}
{
const int boundary_id = 2;
- std::vector<bool> components(n_components, false);
- components[1] = true;
-
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ fe.component_mask(y_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ fe.component_mask(y_displacement));
}
{
const int boundary_id = 4;
- std::vector<bool> components(n_components, false);
- components[2] = true;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ fe.component_mask(z_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ fe.component_mask(z_displacement));
}
{
const int boundary_id = 5;
- std::vector<bool> components(n_components, true);
- components[2] = false;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ (fe.component_mask(x_displacement)
+ |
+ fe.component_mask(y_displacement)));
else
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ (fe.component_mask(x_displacement)
+ |
+ fe.component_mask(y_displacement)));
}
{
const int boundary_id = 6;
- std::vector<bool> components(n_components, true);
- components[2] = false;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ (fe.component_mask(x_displacement)
+ |
+ fe.component_mask(y_displacement)));
else
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
- components);
+ (fe.component_mask(x_displacement)
+ |
+ fe.component_mask(y_displacement)));
}
constraints.close();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
- std::vector<bool> velocity_mask (dim+1+dim, false);
- for (unsigned int d=0; d<dim; ++d)
- velocity_mask[d] = true;
+ const FEValuesExtractors::Vector velocities(0);
VectorTools::interpolate_boundary_values (dof_handler,
1,
StokesBoundaryValues<dim>(),
constraints,
- velocity_mask);
+ fe_collection.component_mask(velocities));
- std::vector<bool> elasticity_mask (dim+1+dim, false);
- for (unsigned int d=dim+1; d<dim+1+dim; ++d)
- elasticity_mask[d] = true;
+ const FEValuesExtractors::Vector displacements(dim+1);
VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction<dim>(dim+1+dim),
constraints,
- elasticity_mask);
+ fe_collection.component_mask(displacements));
}
// There are more constraints we have to
face_q_collection.push_back (stokes_face_quadrature);
face_q_collection.push_back (elasticity_face_quadrature);
- std::vector<bool> stokes_component_mask (dim+1+dim, false);
- for (unsigned int d=0; d<dim; ++d)
- stokes_component_mask[d] = true;
+ const FEValuesExtractors::Vector velocities(0);
KellyErrorEstimator<dim>::estimate (dof_handler,
face_q_collection,
typename FunctionMap<dim>::type(),
solution,
stokes_estimated_error_per_cell,
- stokes_component_mask);
+ fe_collection.component_mask(velocities));
- std::vector<bool> elasticity_component_mask (dim+1+dim, false);
- for (unsigned int d=0; d<dim; ++d)
- elasticity_component_mask[dim+1+d] = true;
+ const FEValuesExtractors::Vector displacements(dim+1);
KellyErrorEstimator<dim>::estimate (dof_handler,
face_q_collection,
typename FunctionMap<dim>::type(),
solution,
elasticity_estimated_error_per_cell,
- elasticity_component_mask);
+ fe_collection.component_mask(displacements));
// We then normalize error estimates by
// dividing by their norm and scale the