]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Reimplement VectorTools::interpolate_based_on_material_id
authorMatthias Maier <tamiko@43-1.org>
Wed, 27 Sep 2017 21:00:01 +0000 (16:00 -0500)
committerMatthias Maier <tamiko@43-1.org>
Thu, 5 Oct 2017 19:38:27 +0000 (14:38 -0500)
include/deal.II/numerics/vector_tools.templates.h

index ada229938abcb0e096863d1165da87a4dd56c337..e15e17fb70894169ad2680712ecff82ed288bbab 100644 (file)
@@ -248,8 +248,10 @@ namespace VectorTools
 
 
     // Internal implementation of interpolate that takes a generic functor
-    // function such that function(cell) is of type
-    // Function<spacedim, typename VectorType::value_type>
+    // function such that function(cell).second is of type
+    // Function<spacedim, typename VectorType::value_type>* and
+    // function(cell).first is a boolean indicating whether the current
+    // cell should be processed
     template <int dim,
               int spacedim,
               typename VectorType,
@@ -370,6 +372,10 @@ namespace VectorTools
           if (fe[fe_index].dofs_per_cell == 0)
             continue;
 
+          // Skip processing of the current cell if we're told to do so.
+          if (!function(cell).first)
+            continue;
+
           // Get transformed, generalized support points
           fe_values.reinit(cell);
           const std::vector<Point<spacedim> > &generalized_support_points =
@@ -390,11 +396,11 @@ namespace VectorTools
           dof_values.resize(n_dofs);
 
           // Get all function values:
-          Assert(n_components == function(cell).n_components,
+          Assert(n_components == function(cell).second->n_components,
                  ExcDimensionMismatch(dof_handler.get_fe().n_components(),
-                                      function(cell).n_components));
-          function(cell).vector_value_list(generalized_support_points,
-                                           function_values);
+                                      function(cell).second->n_components));
+          function(cell).second->vector_value_list(generalized_support_points,
+                                                   function_values);
 
           {
             // Before we can average, we have to transform all function values
@@ -500,9 +506,9 @@ namespace VectorTools
     // internal implementation
     const auto function_map = [&function](
                                 const typename DoFHandlerType<dim, spacedim>::active_cell_iterator &)
-                              -> const Function<spacedim, typename VectorType::value_type> &
+                              -> std::pair<bool, const Function<spacedim, typename VectorType::value_type> *>
     {
-      return function;
+      return std::make_pair(true, &function);
     };
 
     interpolate(mapping, dof_handler, function_map, vec, component_mask);
@@ -591,158 +597,34 @@ namespace VectorTools
   void
   interpolate_based_on_material_id
   (const Mapping<dim,spacedim>        &mapping,
-   const DoFHandlerType<dim,spacedim> &dof,
-   const std::map<types::material_id, const Function<spacedim, typename VectorType::value_type> *> &function_map,
-   VectorType                         &dst,
+   const DoFHandlerType<dim,spacedim> &dof_handler,
+   const std::map<types::material_id, const Function<spacedim, typename VectorType::value_type> *> &functions,
+   VectorType                         &vec,
    const ComponentMask                &component_mask)
   {
-    typedef typename VectorType::value_type number;
-
-    Assert( component_mask.represents_n_components(dof.get_fe(0).n_components()),
-            ExcMessage("The number of components in the mask has to be either "
-                       "zero or equal to the number of components in the finite "
-                       "element.") );
-
-    if ( function_map.size() == 0 )
-      return;
-
-    Assert( function_map.find(numbers::invalid_material_id) == function_map.end(),
-            ExcMessage("You cannot specify the invalid material indicator "
-                       "in your function map."));
-
-    for (typename std::map<types::material_id, const Function<spacedim, number>* >
-         ::const_iterator
-         iter  = function_map.begin();
-         iter != function_map.end();
-         ++iter )
-      {
-        Assert( dof.get_fe(0).n_components() == iter->second->n_components,
-                ExcDimensionMismatch(dof.get_fe(0).n_components(), iter->second->n_components) );
-      }
-
-    const hp::FECollection<dim, spacedim> &fe = dof.get_fe_collection();
-    const unsigned int n_components =  fe.n_components();
-    const bool         fe_is_system = (n_components != 1);
-
-    typename DoFHandlerType<dim,spacedim>::active_cell_iterator
-    cell = dof.begin_active(),
-    endc = dof.end();
-
-    std::vector< std::vector< Point<dim> > > unit_support_points(fe.size());
-    for (unsigned int fe_index = 0; fe_index < fe.size(); ++fe_index)
-      {
-        unit_support_points[fe_index] = fe[fe_index].get_unit_support_points();
-        Assert( unit_support_points[fe_index].size() != 0,
-                ExcNonInterpolatingFE() );
-      }
-
-    std::vector< std::vector<unsigned int> > dofs_of_rep_points(fe.size());
-    std::vector< std::vector<unsigned int> > dof_to_rep_index_table(fe.size());
-    std::vector<unsigned int>                n_rep_points(fe.size(), 0);
-
-    for (unsigned int fe_index = 0; fe_index < fe.size(); ++fe_index)
-      {
-        for (unsigned int i = 0; i < fe[fe_index].dofs_per_cell; ++i)
-          {
-            bool representative = true;
-
-            for (unsigned int j = dofs_of_rep_points[fe_index].size(); j > 0; --j)
-              if ( unit_support_points[fe_index][i]
-                   == unit_support_points[fe_index][dofs_of_rep_points[fe_index][j-1]] )
-                {
-                  dof_to_rep_index_table[fe_index].push_back(j-1);
-                  representative = false;
-                  break;
-                }
-
-            if (representative)
-              {
-                dof_to_rep_index_table[fe_index].push_back(dofs_of_rep_points[fe_index].size());
-                dofs_of_rep_points[fe_index].push_back(i);
-                ++n_rep_points[fe_index];
-              }
-          }
-
-        Assert( dofs_of_rep_points[fe_index].size() == n_rep_points[fe_index],
-                ExcInternalError() );
-        Assert( dof_to_rep_index_table[fe_index].size() == fe[fe_index].dofs_per_cell,
-                ExcInternalError() );
-      }
-
-    const unsigned int max_rep_points = *std::max_element(n_rep_points.begin(),
-                                                          n_rep_points.end());
-    std::vector< types::global_dof_index> dofs_on_cell(fe.max_dofs_per_cell());
-    std::vector< Point<spacedim> > rep_points(max_rep_points);
-
-    std::vector< std::vector<number> >           function_values_scalar(fe.size());
-    std::vector< std::vector< Vector<number> > > function_values_system(fe.size());
-
-    hp::QCollection<dim> support_quadrature;
-    for (unsigned int fe_index = 0; fe_index < fe.size(); ++fe_index)
-      support_quadrature.push_back( Quadrature<dim>(unit_support_points[fe_index]) );
-
-    hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
-    hp::FEValues<dim, spacedim> fe_values(mapping_collection,
-                                          fe,
-                                          support_quadrature,
-                                          update_quadrature_points);
-
-    for ( ; cell != endc; ++cell)
-      if ( cell->is_locally_owned() )
-        if ( function_map.find(cell->material_id()) != function_map.end() )
-          {
-            const unsigned int fe_index = cell->active_fe_index();
-
-            fe_values.reinit(cell);
-
-            const std::vector< Point<spacedim> > &support_points
-              = fe_values.get_present_fe_values().get_quadrature_points();
-
-            rep_points.resize( dofs_of_rep_points[fe_index].size() );
-            for (unsigned int i = 0; i < dofs_of_rep_points[fe_index].size(); ++i)
-              rep_points[i] = support_points[dofs_of_rep_points[fe_index][i]];
-
-            dofs_on_cell.resize( fe[fe_index].dofs_per_cell );
-            cell->get_dof_indices(dofs_on_cell);
-
-            if (fe_is_system)
-              {
-                function_values_system[fe_index].resize( n_rep_points[fe_index],
-                                                         Vector<number>(fe[fe_index].n_components()) );
-
-                function_map.find(cell->material_id())->second->vector_value_list(rep_points,
-                    function_values_system[fe_index]);
-
-                for (unsigned int i = 0; i < fe[fe_index].dofs_per_cell; ++i)
-                  {
-                    const unsigned int component = fe[fe_index].system_to_component_index(i).first;
+    Assert(component_mask.represents_n_components(dof_handler.get_fe().n_components()),
+           ExcMessage("The number of components in the mask has to be either "
+                      "zero or equal to the number of components in the finite "
+                      "element."));
 
-                    if ( component_mask[component] )
-                      {
-                        const unsigned int rep_dof = dof_to_rep_index_table[fe_index][i];
-                        ::dealii::internal::ElementAccess<VectorType>::set(
-                          function_values_system[fe_index][rep_dof](component),
-                          dofs_on_cell[i], dst);
-                      }
-                  }
-              }
-            else
-              {
-                function_values_scalar[fe_index].resize(n_rep_points[fe_index]);
+    Assert (vec.size() == dof_handler.n_dofs(),
+            ExcDimensionMismatch (vec.size(), dof_handler.n_dofs()));
 
-                function_map.find(cell->material_id())->second->value_list
-                (rep_points,
-                 function_values_scalar[fe_index],
-                 0);
+    Assert (component_mask.n_selected_components(dof_handler.get_fe().n_components()) > 0,
+            ComponentMask::ExcNoComponentSelected());
 
-                for (unsigned int i = 0; i < fe[fe_index].dofs_per_cell; ++i)
-                  ::dealii::internal::ElementAccess<VectorType>::set(
-                    function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]],
-                    dofs_on_cell[i], dst);
-              }
-          }
+    // Create a small lambda capture wrapping the function map and call the
+    // internal implementation
+    const auto function_map = [&functions](
+                                const typename DoFHandlerType<dim, spacedim>::active_cell_iterator &cell)
+                              -> std::pair<bool, const Function<spacedim, typename VectorType::value_type> *>
+    {
+      const auto function = functions.find(cell->material_id());
+      bool update = function != functions.end();
+      return std::make_pair(update, function->second);
+    };
 
-    dst.compress (VectorOperation::insert);
+    interpolate(mapping, dof_handler, function_map, vec, component_mask);
   }
 
 

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.