]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Valentin Zingan: Add VectorTools::interpolate_based_on_material_id.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 26 Dec 2013 13:41:56 +0000 (13:41 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 26 Dec 2013 13:41:56 +0000 (13:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@32123 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/numerics/vector_tools.h
deal.II/include/deal.II/numerics/vector_tools.templates.h
deal.II/source/numerics/vector_tools_interpolate.inst.in
tests/deal.II/interpolate_based_on_material_id_01.cc [new file with mode: 0644]
tests/deal.II/interpolate_based_on_material_id_01.output [new file with mode: 0644]

index e2c2969b4b329b582d1e670752ee1f0f4bbd05fb..ccbfba6bc884429617cb8d5ed6018e4018084360 100644 (file)
@@ -85,6 +85,13 @@ inconvenience this causes.
 
 <ol>
 
+  <li> New: The new function VectorTools::interpolate_based_on_material_id()
+  can be used to interpolate several functions onto a mesh, based on the
+  material id of each cell individually.
+  <br>
+  (Valentin Zingan, 2013/12/24)
+  </li>
+
   <li> New: A new reinit() method has been introduced to
   TrilinosWrappers::SparsityPattern that takes all rows that are possibly
   written into as an optional argument. This allows for pre-allocating all
index 9d14cdcc10508727b242ef0805971a85d6e710ab..9558de7a4283858e904b4c0936310b48cb655647 100644 (file)
@@ -502,6 +502,73 @@ namespace VectorTools
                     const InVector           &data_1,
                     OutVector                &data_2);
 
+  /**
+   * This function is a kind of generalization or modification
+   * of the very first interpolate()
+   * function in the series.
+   * It interpolations a set of functions onto the finite element space
+   * given by the DoFHandler argument where the determination which function
+   * to use is made based on the material id (see @ref GlossMaterialId) of
+   * each cell.
+   *
+   * @param mapping        - The mapping to use to determine the location
+   *                         of support points at which the functions
+   *                         are to be evaluated.
+   * @param dof            - DoFHandler initialized with Triangulation
+   *                                                     and
+   *                                                     FiniteElement
+   *                                                     objects,
+   * @param function_map   - std::map reflecting the correspondence
+   *                         between material ids and functions,
+   * @param dst            - global FE vector at the support points,
+   * @param component_mask - mask of components that shall be interpolated
+   *
+   * @note If a material id of some group of cells
+   * is missed in @p function_map, then @p dst will
+   * not be updated in the respective degrees of freedom
+   * of the output vector
+   * For example, if @p dst was successfully
+   * initialized to capture the degrees of freedom of the @p dof_handler
+   * of the problem with all zeros in it,
+   * then those zeros which correspond to
+   * the missed material ids will still remain
+   * in @p dst even after calling this function.
+   *
+   * @note Degrees of freedom located on faces between cells of different
+   * material ids will get their value by that cell which was called
+   * last in the respective loop over cells implemented
+   * in this function.
+   * Since this process is kind of arbitrary,
+   * you cannot control it.
+   * However, if you want to have control over the order in which cells are visited,
+   * let us take a look at the following example: Let @p u be a variable of interest
+   * which is approximated by some CG finite element.
+   * Let @p 0, @p 1 and @p 2 be material ids
+   * of cells on the triangulation.
+   * Let 0: 0.0, 1: 1.0, 2: 2.0
+   * be the whole @p function_map that you want to pass to
+   * this function, where @p key is a material id and
+   * @p value is a value of @p u.
+   * By using the whole @p function_map you do not really know
+   * which values will be assigned to the face DoFs.
+   * On the other hand, if you split the whole @p function_map
+   * into three smaller independent objects
+   * 0: 0.0 and 1: 1.0 and 2: 2.0
+   * and make three distinct calls of this function passing each
+   * of these objects separately (the order depends on what you want
+   * to get between cells), then each subsequent call will rewrite
+   * the intercell @p dofs of the previous one.
+   *
+   * @author Valentin Zingan, 2013
+   */
+  template<typename VECTOR, typename DH>
+  void
+  interpolate_based_on_material_id(const Mapping<DH::dimension, DH::space_dimension>&                          mapping,
+                                   const DH&                                                                   dof_handler,
+                                   const std::map< types::material_id, const Function<DH::space_dimension>* >& function_map,
+                                   VECTOR&                                                                     dst,
+                                   const ComponentMask&                                                        component_mask = ComponentMask());
+
   /**
    * Gives the interpolation of a
    * @p dof1-function @p u1 to a
@@ -2353,6 +2420,12 @@ namespace VectorTools
                              const unsigned int              component);
   //@}
 
+
+  /**
+   * Exception.
+   */
+  DeclException0(ExcInvalidMaterialIndicator);
+
   /**
    * Exception
    */
index 757c2a35fe1eabd90bcba54736478d836ae13b95..cad23ec6b06653216a3251177d818971f3984a6d 100644 (file)
@@ -331,6 +331,154 @@ namespace VectorTools
   }
 
 
+  template<typename VECTOR, typename DH>
+  void
+  interpolate_based_on_material_id(const Mapping<DH::dimension, DH::space_dimension>&                          mapping,
+                                   const DH&                                                                   dof,
+                                   const std::map< types::material_id, const Function<DH::space_dimension>* >& function_map,
+                                   VECTOR&                                                                     dst,
+                                   const ComponentMask&                                                        component_mask)
+  {
+    const unsigned int dim = DH::dimension;
+
+    Assert( component_mask.represents_n_components(dof.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( function_map.size() == 0 )
+      return;
+
+    Assert( function_map.find(numbers::invalid_material_id) == function_map.end(),
+            ExcInvalidMaterialIndicator() );
+
+    for( typename std::map< types::material_id, const Function<DH::space_dimension>* >::const_iterator
+        iter  = function_map.begin();
+        iter != function_map.end();
+        ++iter )
+      {
+        Assert( dof.get_fe().n_components() == iter->second->n_components,
+            ExcDimensionMismatch(dof.get_fe().n_components(), iter->second->n_components) );
+      }
+
+    const hp::FECollection<DH::dimension, DH::space_dimension> fe(dof.get_fe());
+    const unsigned int n_components =  fe.n_components();
+    const bool         fe_is_system = (n_components != 1);
+
+    typename DH::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<unsigned int>                 dofs_on_cell(fe.max_dofs_per_cell());
+    std::vector< Point<DH::space_dimension> > rep_points(max_rep_points);
+
+    std::vector< std::vector<double> >           function_values_scalar(fe.size());
+    std::vector< std::vector< Vector<double> > > 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, DH::space_dimension> mapping_collection(mapping);
+    hp::FEValues<dim, DH::space_dimension> 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<DH::space_dimension> >& 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<double>(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;
+
+                    if( component_mask[component] )
+                      {
+                        const unsigned int rep_dof = dof_to_rep_index_table[fe_index][i];
+                        dst(dofs_on_cell[i])       = function_values_system[fe_index][rep_dof](component);
+                      }
+                  }
+              }
+            else
+              {
+                function_values_scalar[fe_index].resize(n_rep_points[fe_index]);
+
+                function_map.find(cell->material_id())->second->value_list(rep_points,
+                    function_values_scalar[fe_index],
+                    0);
+
+                for(unsigned int i = 0; i < fe[fe_index].dofs_per_cell; ++i)
+                  dst(dofs_on_cell[i]) = function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]];
+              }
+          }
+
+    dst.compress (VectorOperation::insert);
+  }
+
+
   namespace internal
   {
     /**
index bbc0d6b237bbe2ce3ed1ffdfbe33c8ee07a2372e..988939367714de023267dedc809b7acf784ddd86 100644 (file)
@@ -65,6 +65,21 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens
          const FullMatrix<double>&,
          const VEC&,
          VEC&);
+         
+     template
+      void interpolate_based_on_material_id(const Mapping<deal_II_dimension, deal_II_space_dimension>&,
+                                            const DoFHandler<deal_II_dimension, deal_II_space_dimension>&,
+                                            const std::map< types::material_id, const Function<deal_II_space_dimension>* >&,
+                                            VEC&,
+                                            const ComponentMask&);
+         
+     template
+      void interpolate_based_on_material_id(const Mapping<deal_II_dimension>&,
+                                            const hp::DoFHandler<deal_II_dimension>&,
+                                            const std::map< types::material_id, const Function<deal_II_dimension>* >&,
+                                            VEC&,
+                                            const ComponentMask&);
+         
       \}
 #endif
   }
diff --git a/tests/deal.II/interpolate_based_on_material_id_01.cc b/tests/deal.II/interpolate_based_on_material_id_01.cc
new file mode 100644 (file)
index 0000000..64448e1
--- /dev/null
@@ -0,0 +1,123 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2006 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check VectorTools::interpolate_based_on_material_id
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/fe/fe_dgq.h>
+
+#include <fstream>
+#include <vector>
+
+
+template <int dim>
+class F :  public Function<dim>
+{
+public:
+  F (const unsigned int q) : q(q) {}
+
+  virtual double value (const Point<dim> &p,
+                        const unsigned int) const
+  {
+    double v=0;
+    for (unsigned int d=0; d<dim; ++d)
+      for (unsigned int i=0; i<=q; ++i)
+        v += (d+1)*(i+1)*std::pow (p[d], 1.*i);
+    return v;
+  }
+
+private:
+  const unsigned int q;
+};
+
+
+
+template <int dim>
+void test ()
+{
+  Triangulation<dim>     triangulation;
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (3);
+  std::map<types::material_id, const Function<dim>*> functions;
+  for (typename Triangulation<dim>::active_cell_iterator
+        cell = triangulation.begin_active();
+       cell != triangulation.end();
+       ++cell)
+    {
+      cell->set_material_id(cell->index() % 128);
+      if (functions.find(cell->index() % 128) == functions.end())
+       functions[cell->index() % 128]
+         = new ConstantFunction<dim>(cell->index() % 128);
+    }
+  
+  for (unsigned int p=1; p<7-dim; ++p)
+    {
+      FE_DGQ<dim>              fe(p);
+      DoFHandler<dim>        dof_handler(triangulation);
+      dof_handler.distribute_dofs (fe);
+
+      Vector<double> interpolant (dof_handler.n_dofs());
+      VectorTools::interpolate_based_on_material_id (MappingQ1<dim>(),
+                                                    dof_handler,
+                                                    functions,
+                                                    interpolant);
+      for (typename DoFHandler<dim>::active_cell_iterator
+            cell = dof_handler.begin_active();
+          cell != dof_handler.end();
+          ++cell)
+       {
+         Vector<double> values (fe.dofs_per_cell);
+         cell->get_dof_values (interpolant, values);
+         for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+           Assert (values[i] == cell->index() % 128, ExcInternalError());
+       }
+    }
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision (3);
+
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
+
diff --git a/tests/deal.II/interpolate_based_on_material_id_01.output b/tests/deal.II/interpolate_based_on_material_id_01.output
new file mode 100644 (file)
index 0000000..fb71de2
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK

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.