]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert several tutorial programs to the ComponentMask approach.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 27 Sep 2012 16:52:24 +0000 (16:52 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 27 Sep 2012 16:52:24 +0000 (16:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@26791 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-11/step-11.cc
deal.II/examples/step-15/step-15.cc
deal.II/examples/step-17/step-17.cc
deal.II/examples/step-18/step-18.cc
deal.II/examples/step-22/step-22.cc
deal.II/examples/step-31/step-31.cc
deal.II/examples/step-32/step-32.cc
deal.II/examples/step-42/step-42.cc
deal.II/examples/step-43/step-43.cc
deal.II/examples/step-44/step-44.cc
deal.II/examples/step-46/step-46.cc

index 627fcdbbcfecdbe54d532e09be9fabb8383fccef..649cfb1fe0311fca5459430c69b06c89df7d4970 100644 (file)
@@ -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 <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
index c1bf2328bced34b1eed4e1a87dfd844f9c89b1d4..d149599e2bf9d9c7756a07ad4724189a34754229 100644 (file)
@@ -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<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)
index ae10f61ebd4b027c1feb4a34e8d0b894f0a65c0d..a93c3849a830f5cf7bca59a9fb2606c7882b0828 100644 (file)
@@ -1062,7 +1062,7 @@ namespace Step17
                                         typename FunctionMap<dim>::type(),
                                         localized_solution,
                                         local_error_per_cell,
-                                        std::vector<bool>(),
+                                        ComponentMask(),
                                         0,
                                         multithread_info.n_default_threads,
                                         this_mpi_process);
index 58607567163d9e7408531e22e3db47802a53de66..bc051f79c6dea0dbd1c69d1a436d29cb08d47ead 100644 (file)
@@ -1681,8 +1681,7 @@ namespace Step18
                                     // 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,
@@ -1695,7 +1694,7 @@ namespace Step18
                                   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);
@@ -2250,7 +2249,7 @@ namespace Step18
                                        typename FunctionMap<dim>::type(),
                                        incremental_displacement,
                                        error_per_cell,
-                                       std::vector<bool>(),
+                                       ComponentMask(),
                                        0,
                                        multithread_info.n_default_threads,
                                        this_mpi_process);
index 5829d34fbfe9f25cb71c54c883437d2530af2939..20bf102930492a03b43865eb430d6a91006bd355 100644 (file)
@@ -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 <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 &mdash; we
-                                     // first compute the hanging node
+                                     // generated from the DoF handler. Note
+                                     // the order of the two functions &mdash;
+                                     // 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 ();
@@ -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 <int dim>
   void
@@ -1207,14 +1215,13 @@ namespace Step22
   {
     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,
index 0d78aa904670447b53bc0ddd016f12d85d011032..1db324114dcb54bf63f351bbc420428e14fcd24e 100644 (file)
@@ -1469,9 +1469,9 @@ namespace Step31
                          (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;
index e01eb034f470f126f4fe2d232fbafd1eed9ae5a6..f75422b4d8ae2a9b44429c59938e61ff37d28d16 100644 (file)
@@ -2546,13 +2546,12 @@ namespace Step32
       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);
@@ -2858,9 +2857,9 @@ namespace Step32
     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());
@@ -4151,7 +4150,7 @@ namespace Step32
                                         typename FunctionMap<dim>::type(),
                                         temperature_solution,
                                         estimated_error_per_cell,
-                                        std::vector<bool>(),
+                                        ComponentMask(),
                                         0,
                                         0,
                                         triangulation.locally_owned_subdomain());
index 3085e6a0cecb498619bbb646e1d88de0da848adf..1370581865d86d71de7d4ad9bf970c2d1791d995 100644 (file)
@@ -1012,6 +1012,8 @@ namespace Step42
     constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior);
   }
 
+
+
   template <int dim>
   void PlasticityContactProblem<dim>::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<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 ()
   {
@@ -1109,6 +1114,8 @@ namespace Step42
     solution = distributed_solution;
   }
 
+
+
   template <int dim>
   void PlasticityContactProblem<dim>::solve_newton ()
   {
@@ -1119,9 +1126,8 @@ namespace Step42
     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;
@@ -1245,6 +1251,8 @@ namespace Step42
     pcout<< "%%%%%% Rechenzeit output = " << run_time[4] <<std::endl;
   }
 
+
+
   template <int dim>
   void PlasticityContactProblem<dim>::refine_grid ()
   {
@@ -1262,6 +1270,8 @@ namespace Step42
 
   }
 
+
+
   template <int dim>
   void PlasticityContactProblem<dim>::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const
   {
@@ -1293,6 +1303,8 @@ namespace Step42
           }
   }
 
+
+
   template <int dim>
   void PlasticityContactProblem<dim>::output_results (const std::string& title) const
   {
@@ -1351,11 +1363,13 @@ namespace Step42
     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);
 
index 551453329df16cada7be95611f5d38f36cd42773..68763724253c79900110b06a8d2c628fd84b5483 100644 (file)
@@ -821,12 +821,11 @@ namespace Step43
     {
       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 ();
     }
index a60ee65ff43afbe494385c0706e492587fff6a5c..43fc6e54d2979007b20f7bf78714cdc27c76486a 100644 (file)
@@ -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<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();
index 54a73832795d42707f0f6166de3808cbfd14603a..7ecc7c5e84cc9048c223c58a5e43542e75474d41 100644 (file)
@@ -428,23 +428,19 @@ namespace Step46
       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
@@ -1086,25 +1082,21 @@ namespace Step46
     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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.