From df1fd0b1a0ad08e9fce16f65fbcc691537c5312b Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 26 Dec 2013 13:41:56 +0000 Subject: [PATCH] Patch by Valentin Zingan: Add VectorTools::interpolate_based_on_material_id. git-svn-id: https://svn.dealii.org/trunk@32123 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 + .../include/deal.II/numerics/vector_tools.h | 73 +++++++++ .../deal.II/numerics/vector_tools.templates.h | 148 ++++++++++++++++++ .../numerics/vector_tools_interpolate.inst.in | 15 ++ .../interpolate_based_on_material_id_01.cc | 123 +++++++++++++++ ...interpolate_based_on_material_id_01.output | 4 + 6 files changed, 370 insertions(+) create mode 100644 tests/deal.II/interpolate_based_on_material_id_01.cc create mode 100644 tests/deal.II/interpolate_based_on_material_id_01.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index e2c2969b4b..ccbfba6bc8 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -85,6 +85,13 @@ inconvenience this causes.
    +
  1. 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. +
    + (Valentin Zingan, 2013/12/24) +
    2. +
  2. 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 diff --git a/deal.II/include/deal.II/numerics/vector_tools.h b/deal.II/include/deal.II/numerics/vector_tools.h index 9d14cdcc10..9558de7a42 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.h +++ b/deal.II/include/deal.II/numerics/vector_tools.h @@ -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 + void + interpolate_based_on_material_id(const Mapping& mapping, + const DH& dof_handler, + const std::map< types::material_id, const Function* >& 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 */ diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 757c2a35fe..cad23ec6b0 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -331,6 +331,154 @@ namespace VectorTools } + template + void + interpolate_based_on_material_id(const Mapping& mapping, + const DH& dof, + const std::map< types::material_id, const Function* >& 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* >::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 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 > > 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 > dofs_of_rep_points(fe.size()); + std::vector< std::vector > dof_to_rep_index_table(fe.size()); + std::vector 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 dofs_on_cell(fe.max_dofs_per_cell()); + std::vector< Point > rep_points(max_rep_points); + + std::vector< std::vector > function_values_scalar(fe.size()); + std::vector< std::vector< Vector > > function_values_system(fe.size()); + + hp::QCollection support_quadrature; + for(unsigned int fe_index = 0; fe_index < fe.size(); ++fe_index) + support_quadrature.push_back( Quadrature(unit_support_points[fe_index]) ); + + hp::MappingCollection mapping_collection(mapping); + hp::FEValues 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 >& 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(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 { /** diff --git a/deal.II/source/numerics/vector_tools_interpolate.inst.in b/deal.II/source/numerics/vector_tools_interpolate.inst.in index bbc0d6b237..9889393677 100644 --- a/deal.II/source/numerics/vector_tools_interpolate.inst.in +++ b/deal.II/source/numerics/vector_tools_interpolate.inst.in @@ -65,6 +65,21 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const FullMatrix&, const VEC&, VEC&); + + template + void interpolate_based_on_material_id(const Mapping&, + const DoFHandler&, + const std::map< types::material_id, const Function* >&, + VEC&, + const ComponentMask&); + + template + void interpolate_based_on_material_id(const Mapping&, + const hp::DoFHandler&, + const std::map< types::material_id, const Function* >&, + 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 index 0000000000..64448e1d2a --- /dev/null +++ b/tests/deal.II/interpolate_based_on_material_id_01.cc @@ -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 +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +class F : public Function +{ +public: + F (const unsigned int q) : q(q) {} + + virtual double value (const Point &p, + const unsigned int) const + { + double v=0; + for (unsigned int d=0; d +void test () +{ + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (3); + std::map*> functions; + for (typename Triangulation::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(cell->index() % 128); + } + + for (unsigned int p=1; p<7-dim; ++p) + { + FE_DGQ fe(p); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs (fe); + + Vector interpolant (dof_handler.n_dofs()); + VectorTools::interpolate_based_on_material_id (MappingQ1(), + dof_handler, + functions, + interpolant); + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); + ++cell) + { + Vector values (fe.dofs_per_cell); + cell->get_dof_values (interpolant, values); + for (unsigned int i=0; iindex() % 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 index 0000000000..fb71de2867 --- /dev/null +++ b/tests/deal.II/interpolate_based_on_material_id_01.output @@ -0,0 +1,4 @@ + +DEAL::OK +DEAL::OK +DEAL::OK -- 2.39.5