<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
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
const unsigned int component);
//@}
+
+ /**
+ * Exception.
+ */
+ DeclException0(ExcInvalidMaterialIndicator);
+
/**
* Exception
*/
}
+ 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
{
/**
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
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $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>();
+}
+
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
projects / dealii-svn.git / commitdiff