]> https://gitweb.dealii.org/ - dealii.git/commitdiff
intermediate
authorMatthias Maier <tamiko@43-1.org>
Sat, 2 Sep 2017 04:42:07 +0000 (23:42 -0500)
committerMatthias Maier <tamiko@43-1.org>
Sun, 17 Sep 2017 20:26:40 +0000 (15:26 -0500)
include/deal.II/fe/fe.h
include/deal.II/numerics/vector_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/numerics/vector_tools_interpolate.inst.in

index 62f0ae008517111b3148bef2a056eb14adc7697a..b55fdcec1a43bf534dd86cb426ad22a29702f30a 100644 (file)
@@ -1972,6 +1972,7 @@ public:
    */
   bool has_generalized_support_points () const;
 
+//FIXME
   /**
    * Return the equivalent to get_generalized_support_points(), except
    * for faces.
@@ -1979,6 +1980,7 @@ public:
   const std::vector<Point<dim-1> > &
   get_generalized_face_support_points () const;
 
+//FIXME
   /**
    * Return whether a finite element has defined generalized support points on
    * faces. If the result is true, then a call to the
index b7658cb197295c7617d2042054920c778fb00d0a..cd3a3536e9e55b71ba80f61be02f5f0f9fa2f4fb 100644 (file)
@@ -566,6 +566,9 @@ namespace VectorTools
    * @name Interpolation and projection
    */
   //@{
+
+
+
   /**
    * Compute the interpolation of @p function at the support points to the
    * finite element space described by the Triangulation and FiniteElement
index 9ff7b90ff5e8fe86808efa560e3d67b04a714602..27adfdd4fd60fdc9748a02dac778c7bd63fbf455 100644 (file)
@@ -85,215 +85,282 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace VectorTools
 {
-  template <int dim, int spacedim, typename VectorType,
+  template <int dim,
+            int spacedim,
+            typename VectorType,
             template <int, int> class DoFHandlerType>
-  void interpolate (const Mapping<dim,spacedim>                               &mapping,
-                    const DoFHandlerType<dim,spacedim>                        &dof,
-                    const Function<spacedim, typename VectorType::value_type> &function,
-                    VectorType                                                &vec,
-                    const ComponentMask                                       &component_mask)
+  void interpolate(
+    const Mapping<dim, spacedim> &mapping,
+    const DoFHandlerType<dim, spacedim> &dof_handler,
+    const Function<spacedim, typename VectorType::value_type> &function,
+    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.") );
+    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."));
 
-    Assert (vec.size() == dof.n_dofs(),
-            ExcDimensionMismatch (vec.size(), dof.n_dofs()));
-    Assert (dof.get_fe(0).n_components() == function.n_components,
-            ExcDimensionMismatch(dof.get_fe(0).n_components(),
-                                 function.n_components));
-    Assert (component_mask.n_selected_components(dof.get_fe(0).n_components()) > 0,
-            ComponentMask::ExcNoComponentSelected());
+    Assert (vec.size() == dof_handler.n_dofs(),
+            ExcDimensionMismatch (vec.size(), dof_handler.n_dofs()));
 
-    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);
+    Assert (dof_handler.get_fe().n_components() == function.n_components,
+            ExcDimensionMismatch(dof_handler.get_fe().n_components(), function.n_components));
 
-    typename DoFHandlerType<dim,spacedim>::active_cell_iterator
-    cell = dof.begin_active(),
-    endc = dof.end();
+    Assert (component_mask.n_selected_components(dof_handler.get_fe().n_components()) > 0,
+            ComponentMask::ExcNoComponentSelected());
 
-    // For FESystems many of the
-    // unit_support_points will appear
-    // multiple times, as a point may be
-    // unit_support_point for several of the
-    // components of the system.  The following
-    // is rather complicated, but at least
-    // attempts to avoid evaluating the
-    // vectorfunction multiple times at the
-    // same point on a cell.
     //
-    // First check that the desired components are interpolating.
+    // Computing the generalized interpolant isn't quite as straightforward
+    // as for classical Lagrange elements. A major complication is the fact
+    // it generally doesn't hold true that a function evaluates to the same
+    // dof coefficient on different cells. This means *setting* the value
+    // of a (global) degree of freedom computed on one cell doesn't
+    // necessarily lead to the same result when computed on a neighboring
+    // cell (that shares the same global degree of freedom).
+    //
+    // We thus, do the following operation:
+    //
+    // On each cell:
+    //
+    //  - We first determine all function values u(x_i) in generalized
+    //    support points
+    //
+    //  - We transform these function values back to the unit cell
+    //    according to the conformity of the component (scalar, Hdiv, or
+    //    Hcurl conforming); see [Monk, Finite Element Methods for Maxwell's
+    //    Equations, p.77ff Section 3.9] for details. This results in
+    //    \hat u(\hat x_i)
+    //
+    //  - We convert these generalized support point values to nodal values
+    //
+    //  - For every global dof we take the average 1 / n_K \sum_{K} dof_K
+    //    where n_K is the number of cells sharing the global dof and dof_K
+    //    is the computed value on the cell K.
+    //
+    // For every degree of freedom that is shared by k cells, we compute
+    // its value on all k cells and take the weighted average with respect
+    // to the JxW values.
+    //
+
+    typedef typename VectorType::value_type number;
+
+    const hp::FECollection<dim, spacedim> fe(dof_handler.get_fe());
+
+    std::vector<types::global_dof_index> dofs_on_cell(fe.max_dofs_per_cell());
+
+    // Temporary storage for cell-wise interpolation operation. We store a
+    // variant for every fe we encounter to speed up resizing operations.
+    // The first vector is used for local function evaluation. The vector
+    // dof_values is used to store intermediate cell-wise interpolation
+    // results (see the detailed explanation in the for loop further down
+    // below).
+
+    std::vector<std::vector<Vector<number> > > fe_function_values(fe.size());
+    std::vector<std::vector<number> > fe_dof_values(fe.size());
+
+    // We will need two temporary global vectors that store the new values
+    // and weights.
+    VectorType interpolation;
+    VectorType weights;
+    interpolation.reinit(vec);
+    weights.reinit(vec);
+
+    // We use an FEValues object to transform all generalized support
+    // points from the unit cell to the real cell coordinates. Thus,
+    // initialize a quadrature with all generalized support points and
+    // create an FEValues object with it.
+
+    hp::QCollection<dim> support_quadrature;
     for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
       {
-        for (unsigned int component_index = 0; component_index < n_components; ++component_index)
-          {
-            if (component_mask[component_index] == true)
-              {
-                Assert ((fe[fe_index].base_element(fe[fe_index].component_to_base_index
-                                                   (component_index).first).has_support_points()) ||
-                        (fe[fe_index].base_element(fe[fe_index].component_to_base_index
-                                                   (component_index).first).dofs_per_cell == 0),
-                        ExcNonInterpolatingFE());
-              }
-          }
+        const auto &points = fe[fe_index].get_generalized_support_points();
+        support_quadrature.push_back(Quadrature<dim>(points));
       }
 
-    // Find the support points on a cell that are mentioned multiple times, and
-    // only add each once.  Each multiple point gets to know the dof index of
-    // its representative point by the dof_to_rep_dof_table.
+    hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
 
-    // the following vector collects all unit support points p[i],
-    // 0<=i<fe.dofs_per_cell, for which unit_support_point(i) is unique
-    // (representative). The position of a support point within this vector is
-    // called the rep index.
-    std::vector<std::vector<Point<dim> > > rep_unit_support_points (fe.size());
-    // the following table converts a dof i
-    // to the rep index.
-    std::vector<std::vector<types::global_dof_index> > dof_to_rep_index_table(fe.size());
+    hp::FEValues<dim, spacedim> fe_values(
+      mapping_collection,
+      fe,
+      support_quadrature,
+      update_quadrature_points);
+
+    // An extra FEValues object to compute jacobians.
+    // Only re-initialized in case of Hcurl or Hdiv conforming elements,
+    // i.e. if we really need the information.
+    hp::FEValues<dim, spacedim> fe_values_jacobians(
+      mapping_collection,
+      fe,
+      support_quadrature,
+      update_jacobians | update_inverse_jacobians);
 
-    std::vector<unsigned int> n_rep_points (fe.size(), 0);
+    //
+    // Now loop over all locally owned, active cells.
+    //
 
-    for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
+    for (auto cell : dof_handler.active_cell_iterators())
       {
-        for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
+        // If this cell is not locally owned, do nothing.
+        if (! cell->is_locally_owned())
+          continue;
+
+        const unsigned int fe_index = cell->active_fe_index();
+
+        // Do nothing if there are no local degrees of freedom.
+        if (fe[fe_index].dofs_per_cell == 0)
+          continue;
+
+        // Get transformed, generalized support points
+        fe_values.reinit(cell);
+        const std::vector<Point<spacedim> > &generalized_support_points =
+          fe_values.get_present_fe_values().get_quadrature_points();
+
+        // Get indices of the dofs on this cell
+        const auto n_dofs = fe[fe_index].dofs_per_cell;
+        dofs_on_cell.resize (n_dofs);
+        cell->get_dof_indices (dofs_on_cell);
+
+        // Prepare temporary storage
+        auto &function_values = fe_function_values[fe_index];
+        auto &dof_values = fe_dof_values[fe_index];
+
+        const auto n_components = fe[fe_index].n_components();
+        function_values.resize(generalized_support_points.size(),
+                               Vector<number>(n_components));
+        dof_values.resize(n_dofs);
+
+        // Get all function values:
+        function.vector_value_list(generalized_support_points,
+                                   function_values);
+
+        // FIXME: In case of an FESystem we have to apply this
+        // transformation according to the conformity of each base element.
+
+        // Before we can average, we have to transform all function values
+        // from the real cell back to the unit cell. We query the finite
+        // element for the correct transformation.
+        switch (fe[fe_index].conforming_space)
           {
-            const unsigned int component
-              = fe[fe_index].system_to_component_index(i).first;
-            if (component_mask[component] == true)
+          case FiniteElementData<dim>::H1:
+            DEAL_II_FALLTHROUGH;
+          case FiniteElementData<dim>::L2:
+            // See Monk, Finite Element Methods for Maxwell's Equations,
+            // p. 77ff, formula (3.74).
+            // For given mapping F_K: \hat K \to K, we have to transform
+            //  \hat p = p\circ F_K
+            //  i.e., do nothing.
+            break;
+
+          case FiniteElementData<dim>::Hcurl:
+            // See Monk, Finite Element Methods for Maxwell's Equations,
+            // p. 77ff, formula (3.76) and Corollary 3.58.
+            // For given mapping F_K: \hat K \to K, we have to transform
+            //  \hat u = (dF_K)^T u\circ F_K
+
+            fe_values_jacobians.reinit(cell);
+            for (unsigned int i = 0; i < function_values.size(); ++i)
               {
-                const Point<dim> dof_support_point = fe[fe_index].unit_support_point(i);
-
-                bool representative=true;
-                // the following loop is looped
-                // the other way round to get
-                // the minimal effort of
-                // O(fe.dofs_per_cell) for multiple
-                // support points that are placed
-                // one after the other.
-                for (unsigned int j=rep_unit_support_points[fe_index].size(); j>0; --j)
-                  if (dof_support_point
-                      == rep_unit_support_points[fe_index][j-1])
-                    {
-                      dof_to_rep_index_table[fe_index].push_back(j-1);
-                      representative=false;
-                      break;
-                    }
+                const auto &jacobians =
+                  fe_values_jacobians.get_present_fe_values().get_jacobians();
 
-                if (representative)
-                  {
-                    dof_to_rep_index_table[fe_index].push_back
-                    (rep_unit_support_points[fe_index].size());
-                    rep_unit_support_points[fe_index].push_back(dof_support_point);
-                    ++n_rep_points[fe_index];
-                  }
+                // value[m] <- sum jacobian_transpose[m][n] * old_value[n]:
+                const auto old_value = function_values[i];
+                TensorAccessors::contract<1, 2, 1, dim>(
+                  function_values[i], jacobians[i].transpose(), old_value);
               }
-            else
+            break;
+
+          case FiniteElementData<dim>::Hdiv:
+            // See Monk, Finite Element Methods for Maxwell's Equations,
+            // p. 79ff, formula (3.77) and Lemma 3.59.
+            // For given mapping F_K: \hat K \to K, we have to transform
+            //  \hat w = det(dF_K) (dF_K)^{-1} w\circ F_K
+
+            fe_values_jacobians.reinit(cell);
+            for (unsigned int i = 0; i < function_values.size(); ++i)
               {
-                // If correct component not to be interpolated, append invalid index to in table
-                dof_to_rep_index_table[fe_index].push_back(numbers::invalid_dof_index);
+                const auto &jacobians =
+                  fe_values_jacobians.get_present_fe_values().get_jacobians();
+                const auto &inverse_jacobians =
+                  fe_values_jacobians.get_present_fe_values().get_inverse_jacobians();
+
+                // value[m] <- sum inverse_jacobians[m][n] * old_value[n]:
+                const auto old_value = function_values[i];
+                TensorAccessors::contract<1, 2, 1, dim>(
+                  function_values[i], inverse_jacobians[i], old_value);
+
+                for (unsigned int j = 0; j < n_components; ++j)
+                  function_values[i][j] *= jacobians[i].determinant();
               }
-          }
-
-        Assert(rep_unit_support_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());
-      }
+            break;
 
-    std::vector<types::global_dof_index> dofs_on_cell (fe.max_dofs_per_cell());
+          default:
+            Assert(false,
+                   ExcMessage(
+                     "The supplied finite element has an unknown conformity."));
+          } /*switch*/
 
-    // get space for the values of the
-    // function at the rep support points.
-    //
-    // have two versions, one for system fe
-    // and one for scalar ones, to take the
-    // more efficient one respectively
-    std::vector<std::vector<number> >         function_values_scalar(fe.size());
-    std::vector<std::vector<Vector<number> > > function_values_system(fe.size());
-
-    // Make a quadrature rule from support points
-    // to feed it into FEValues
-    hp::QCollection<dim> support_quadrature;
-    for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
-      support_quadrature.push_back (Quadrature<dim>(rep_unit_support_points[fe_index]));
+        FETools::convert_generalized_support_point_values_to_dof_values(
+          fe[fe_index], function_values, dof_values);
 
-    // Transformed support points are computed by
-    // FEValues
-    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
+        for (unsigned int i=0; i < n_dofs; ++i)
+          {
+            ::dealii::internal::ElementAccess<VectorType>::add(
+              typename VectorType::value_type(1.0),
+              dofs_on_cell[i],
+              weights);
 
-    hp::FEValues<dim,spacedim> fe_values (mapping_collection,
-                                          fe, support_quadrature, update_quadrature_points);
+            const auto &nonzero_components =
+              fe[fe_index].get_nonzero_components(i);
 
-    for (; cell!=endc; ++cell)
-      if (cell->is_locally_owned())
-        {
-          const unsigned int fe_index = cell->active_fe_index();
-          if (fe[fe_index].dofs_per_cell != 0)
-            {
-              // for each cell:
-              // get location of finite element
-              // support_points
-              fe_values.reinit(cell);
-              const std::vector<Point<spacedim> > &rep_support_points =
-                fe_values.get_present_fe_values().get_quadrature_points();
+            // Figure out whether the component mask applies. We assume
+            // that we are allowed to set degrees of freedom if at least
+            // one of the components (of the dof) is selected.
+            //
+            // FIXME: Add a debug assert.
+            bool selected = false;
+            for (unsigned int i = 0; i < nonzero_components.size(); ++i)
+              selected =
+                selected || (nonzero_components[i] && component_mask[i]);
 
-              // get indices of the dofs on this cell
-              dofs_on_cell.resize (fe[fe_index].dofs_per_cell);
-              cell->get_dof_indices (dofs_on_cell);
+            if (selected)
+              {
+                // Add local values to the global vectors
+                ::dealii::internal::ElementAccess<VectorType>::add(
+                  dof_values[i], dofs_on_cell[i], interpolation);
+              }
+            else
+              {
+                // If a component is ignored, simply copy all dof values
+                // from the vector "vec":
+                const auto value =
+                  ::dealii::internal::ElementAccess<VectorType>::get(
+                    vec, dofs_on_cell[i]);
+                ::dealii::internal::ElementAccess<VectorType>::add(
+                  value, dofs_on_cell[i], interpolation);
+              }
+          }
+      } /* loop over dof_handler.active_cell_iterators() */
 
+    interpolation.compress(VectorOperation::add);
+    weights.compress(VectorOperation::add);
 
-              if (fe_is_system)
-                {
-                  // get function values at
-                  // these points. Here: get
-                  // all components
-                  function_values_system[fe_index]
-                  .resize (n_rep_points[fe_index],
-                           Vector<number> (fe[fe_index].n_components()));
-                  function.vector_value_list (rep_support_points,
-                                              function_values_system[fe_index]);
-                  // distribute the function
-                  // values to the global
-                  // vector
-                  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;
-                      if (component_mask[component] == true)
-                        {
-                          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], vec);
-                        }
-                    }
-                }
-              else
-                {
-                  // get first component only,
-                  // which is the only component
-                  // in the function anyway
-                  function_values_scalar[fe_index].resize (n_rep_points[fe_index]);
-                  function.value_list (rep_support_points,
-                                       function_values_scalar[fe_index],
-                                       0);
-                  // distribute the function
-                  // values to the global
-                  // vector
-                  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], vec);
-                }
-            }
-        }
+    for (const auto i : interpolation.locally_owned_elements())
+      {
+        const auto value =
+          ::dealii::internal::ElementAccess<VectorType>::get(interpolation, i);
+        const auto weight =
+          ::dealii::internal::ElementAccess<VectorType>::get(weights, i);
+        ::dealii::internal::ElementAccess<VectorType>::set(
+          value / weight, i, vec);
+      }
     vec.compress(VectorOperation::insert);
   }
 
 
+
+
   template <int dim, int spacedim, typename VectorType,
             template <int,int> class DoFHandlerType>
   void interpolate (const DoFHandlerType<dim,spacedim>                       &dof,
index 6e1d4d615a97b229fc2490e31db07e95c1bacd3e..b8e1bf796446e5fef77897336a3bf2dd1f31c8f5 100644 (file)
@@ -14,6 +14,7 @@
 // ---------------------------------------------------------------------
 
 
+
 for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
 {
 #if deal_II_dimension <= deal_II_space_dimension

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.