From: Wolfgang Bangerth Date: Thu, 27 Sep 2012 16:52:24 +0000 (+0000) Subject: Convert several tutorial programs to the ComponentMask approach. X-Git-Tag: v8.0.0~2018 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=578ab67811c1168746f062ac5b1bac616b0f6386;p=dealii.git Convert several tutorial programs to the ComponentMask approach. git-svn-id: https://svn.dealii.org/trunk@26791 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-11/step-11.cc b/deal.II/examples/step-11/step-11.cc index 627fcdbbcf..649cfb1fe0 100644 --- a/deal.II/examples/step-11/step-11.cc +++ b/deal.II/examples/step-11/step-11.cc @@ -173,13 +173,24 @@ namespace Step11 // 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 true. + // anyway, this mask in reality should + // have only one entry with a + // true value. However, + // the ComponentMask class has + // semantics that allow it to + // represents a mask of indefinite + // size whose every element equals + // true when one just + // default constructs such an object, + // so this is what we'll do here. std::vector boundary_dofs (dof_handler.n_dofs(), false); - DoFTools::extract_boundary_dofs (dof_handler, std::vector(1,true), + DoFTools::extract_boundary_dofs (dof_handler, + ComponentMask(), boundary_dofs); // Now first for the generation of diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index c1bf2328bc..d149599e2b 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -702,12 +702,13 @@ namespace Step15 // 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 boundary_dofs (dof_handler.n_dofs()); DoFTools::extract_boundary_dofs (dof_handler, - std::vector(1,true), + ComponentMask(), boundary_dofs); for (unsigned int i=0; i::type(), localized_solution, local_error_per_cell, - std::vector(), + ComponentMask(), 0, multithread_info.n_default_threads, this_mpi_process); diff --git a/deal.II/examples/step-18/step-18.cc b/deal.II/examples/step-18/step-18.cc index 5860756716..bc051f79c6 100644 --- a/deal.II/examples/step-18/step-18.cc +++ b/deal.II/examples/step-18/step-18.cc @@ -1681,8 +1681,7 @@ namespace Step18 // restrict vertical motion, it // has only its last component // set: - std::vector z_component (dim, false); - z_component[dim-1] = true; + FEValuesExtractors::Scalar z_component (dim-1); std::map boundary_values; VectorTools:: interpolate_boundary_values (dof_handler, @@ -1695,7 +1694,7 @@ namespace Step18 IncrementalBoundaryValues(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); @@ -2250,7 +2249,7 @@ namespace Step18 typename FunctionMap::type(), incremental_displacement, error_per_cell, - std::vector(), + ComponentMask(), 0, multithread_info.n_default_threads, this_mpi_process); diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index 5829d34fbf..20bf102930 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -573,52 +573,55 @@ namespace Step22 // 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 component_mask 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 - // H1 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 H1 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 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(), constraints, - component_mask); + fe.component_mask(velocities)); } constraints.close (); @@ -1198,8 +1201,13 @@ namespace Step22 // 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 void @@ -1207,14 +1215,13 @@ namespace Step22 { Vector estimated_error_per_cell (triangulation.n_active_cells()); - std::vector component_mask (dim+1, false); - component_mask[dim] = true; + FEValuesExtractors::Scalar pressure(dim); KellyErrorEstimator::estimate (dof_handler, QGauss(degree+1), typename FunctionMap::type(), solution, estimated_error_per_cell, - component_mask); + fe.component_mask(pressure)); GridRefinement::refine_and_coarsen_fixed_number (triangulation, estimated_error_per_cell, diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 0d78aa9046..1db324114d 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -1469,9 +1469,9 @@ namespace Step31 (new TrilinosWrappers::PreconditionAMG()); std::vector > constant_modes; - std::vector 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; diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index e01eb034f4..f75422b4d8 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -2546,13 +2546,12 @@ namespace Step32 DoFTools::make_hanging_node_constraints (stokes_dof_handler, stokes_constraints); - std::vector velocity_mask (dim+1, true); - velocity_mask[dim] = false; + FEValuesExtractors::Vector velocity_components(0); VectorTools::interpolate_boundary_values (stokes_dof_handler, 0, ZeroFunction(dim+1), stokes_constraints, - velocity_mask); + stokes_fe.component_mask(velocity_components)); std::set no_normal_flux_boundaries; no_normal_flux_boundaries.insert (1); @@ -2858,9 +2857,9 @@ namespace Step32 assemble_stokes_preconditioner (); std::vector > constant_modes; - std::vector 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()); @@ -4151,7 +4150,7 @@ namespace Step32 typename FunctionMap::type(), temperature_solution, estimated_error_per_cell, - std::vector(), + ComponentMask(), 0, 0, triangulation.locally_owned_subdomain()); diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 3085e6a0ce..1370581865 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -1012,6 +1012,8 @@ namespace Step42 constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior); } + + template void PlasticityContactProblem::dirichlet_constraints () { @@ -1028,27 +1030,30 @@ namespace Step42 constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs); constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes); - std::vector 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(), 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(), constraints_dirichlet_hanging_nodes, - component_mask); + (fe.component_mask(x_displacement) + | + fe.component_mask(y_displacement))); constraints_dirichlet_hanging_nodes.close (); } + + template void PlasticityContactProblem::solve () { @@ -1109,6 +1114,8 @@ namespace Step42 solution = distributed_solution; } + + template void PlasticityContactProblem::solve_newton () { @@ -1119,9 +1126,8 @@ namespace Step42 Timer t; std::vector > constant_modes; - std::vector 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; @@ -1245,6 +1251,8 @@ namespace Step42 pcout<< "%%%%%% Rechenzeit output = " << run_time[4] < void PlasticityContactProblem::refine_grid () { @@ -1262,6 +1270,8 @@ namespace Step42 } + + template void PlasticityContactProblem::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const { @@ -1293,6 +1303,8 @@ namespace Step42 } } + + template void PlasticityContactProblem::output_results (const std::string& title) const { @@ -1351,11 +1363,13 @@ namespace Step42 move_mesh (tmp); } + + template void PlasticityContactProblem::run () { pcout << "Solving problem in " << dim << " space dimensions." << std::endl; - + Timer t; run_time.resize (8); diff --git a/deal.II/examples/step-43/step-43.cc b/deal.II/examples/step-43/step-43.cc index 551453329d..6876372425 100644 --- a/deal.II/examples/step-43/step-43.cc +++ b/deal.II/examples/step-43/step-43.cc @@ -821,12 +821,11 @@ namespace Step43 { darcy_preconditioner_constraints.clear (); - std::vector 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 (); } diff --git a/deal.II/examples/step-44/step-44.cc b/deal.II/examples/step-44/step-44.cc index a60ee65ff4..43fc6e54d2 100644 --- a/deal.II/examples/step-44/step-44.cc +++ b/deal.II/examples/step-44/step-44.cc @@ -2641,97 +2641,114 @@ namespace Step44 // +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 components(n_components, false); - components[0] = true; - if (apply_dirichlet_bc == true) VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, - components); + fe.component_mask(x_displacement)); else VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, - components); + fe.component_mask(x_displacement)); } { const int boundary_id = 2; - std::vector components(n_components, false); - components[1] = true; - if (apply_dirichlet_bc == true) VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, - components); + fe.component_mask(y_displacement)); else VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, - components); + fe.component_mask(y_displacement)); } { const int boundary_id = 4; - std::vector components(n_components, false); - components[2] = true; if (apply_dirichlet_bc == true) VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, - components); + fe.component_mask(z_displacement)); else VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, - components); + fe.component_mask(z_displacement)); } { const int boundary_id = 5; - std::vector components(n_components, true); - components[2] = false; if (apply_dirichlet_bc == true) VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(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(n_components), constraints, - components); + (fe.component_mask(x_displacement) + | + fe.component_mask(y_displacement))); } { const int boundary_id = 6; - std::vector components(n_components, true); - components[2] = false; if (apply_dirichlet_bc == true) VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(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(n_components), constraints, - components); + (fe.component_mask(x_displacement) + | + fe.component_mask(y_displacement))); } constraints.close(); diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index 54a7383279..7ecc7c5e84 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -428,23 +428,19 @@ namespace Step46 DoFTools::make_hanging_node_constraints (dof_handler, constraints); - std::vector velocity_mask (dim+1+dim, false); - for (unsigned int d=0; d(), constraints, - velocity_mask); + fe_collection.component_mask(velocities)); - std::vector elasticity_mask (dim+1+dim, false); - for (unsigned int d=dim+1; d(dim+1+dim), constraints, - elasticity_mask); + fe_collection.component_mask(displacements)); } // There are more constraints we have to @@ -1086,25 +1082,21 @@ namespace Step46 face_q_collection.push_back (stokes_face_quadrature); face_q_collection.push_back (elasticity_face_quadrature); - std::vector stokes_component_mask (dim+1+dim, false); - for (unsigned int d=0; d::estimate (dof_handler, face_q_collection, typename FunctionMap::type(), solution, stokes_estimated_error_per_cell, - stokes_component_mask); + fe_collection.component_mask(velocities)); - std::vector elasticity_component_mask (dim+1+dim, false); - for (unsigned int d=0; d::estimate (dof_handler, face_q_collection, typename FunctionMap::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