From fad0dfb7a9412258ee4a2f8910bebbf6c55316ed Mon Sep 17 00:00:00 2001 From: Matthias Maier <tamiko@43-1.org> Date: Mon, 29 Feb 2016 13:41:43 -0600 Subject: [PATCH] take over #2273 from tjhei/split_instantiations Split instantiations of matrix_tools.cc Conflicts: source/numerics/matrix_tools.cc source/numerics/matrix_tools.inst.in --- .../numerics/matrix_creator.templates.h | 1894 ++++++++++++++ source/numerics/CMakeLists.txt | 5 + source/numerics/matrix_creator.cc | 28 + source/numerics/matrix_creator.inst.in | 315 +++ source/numerics/matrix_creator_inst2.cc | 20 + source/numerics/matrix_creator_inst3.cc | 20 + source/numerics/matrix_tools.cc | 2315 ----------------- source/numerics/matrix_tools.inst.in | 309 +-- source/numerics/matrix_tools_once.cc | 520 ++++ 9 files changed, 2815 insertions(+), 2611 deletions(-) create mode 100644 include/deal.II/numerics/matrix_creator.templates.h create mode 100644 source/numerics/matrix_creator.cc create mode 100644 source/numerics/matrix_creator.inst.in create mode 100644 source/numerics/matrix_creator_inst2.cc create mode 100644 source/numerics/matrix_creator_inst3.cc create mode 100644 source/numerics/matrix_tools_once.cc diff --git a/include/deal.II/numerics/matrix_creator.templates.h b/include/deal.II/numerics/matrix_creator.templates.h new file mode 100644 index 0000000000..30953c19e3 --- /dev/null +++ b/include/deal.II/numerics/matrix_creator.templates.h @@ -0,0 +1,1894 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__matrix_creator_templates_h +#define dealii__matrix_creator_templates_h + +#include <deal.II/base/function.h> +#include <deal.II/base/quadrature.h> +#include <deal.II/base/work_stream.h> +#include <deal.II/base/geometry_info.h> +#include <deal.II/base/quadrature.h> +#include <deal.II/dofs/dof_handler.h> +#include <deal.II/dofs/dof_accessor.h> +#include <deal.II/dofs/dof_tools.h> +#include <deal.II/fe/fe.h> +#include <deal.II/fe/fe_values.h> +#include <deal.II/fe/mapping_q1.h> +#include <deal.II/grid/tria_iterator.h> +#include <deal.II/hp/fe_values.h> +#include <deal.II/hp/mapping_collection.h> +#include <deal.II/numerics/matrix_tools.h> +#include <deal.II/lac/vector.h> +#include <deal.II/lac/block_vector.h> +#include <deal.II/lac/full_matrix.h> +#include <deal.II/lac/sparse_matrix.h> +#include <deal.II/lac/block_sparse_matrix.h> + +#ifdef DEAL_II_WITH_PETSC +# include <deal.II/lac/petsc_parallel_sparse_matrix.h> +# include <deal.II/lac/petsc_sparse_matrix.h> +# include <deal.II/lac/petsc_parallel_vector.h> +# include <deal.II/lac/petsc_vector.h> +# include <deal.II/lac/petsc_parallel_block_sparse_matrix.h> +#endif + +#ifdef DEAL_II_WITH_TRILINOS +# include <deal.II/lac/trilinos_sparse_matrix.h> +# include <deal.II/lac/trilinos_vector.h> +# include <deal.II/lac/trilinos_block_sparse_matrix.h> +# include <deal.II/lac/trilinos_block_vector.h> +#endif + +#include <algorithm> + + +#include <algorithm> +#include <set> +#include <cmath> + + +DEAL_II_NAMESPACE_OPEN + +namespace MatrixCreator +{ + namespace internal + { + /** + * Convenience abbreviation for + * pairs of DoF handler cell + * iterators. This type works + * just like a + * <tt>std::pair<iterator,iterator></tt> + * but is templatized on the + * dof handler that should be used. + */ + template <typename DoFHandlerType> + struct IteratorRange + { + /** + * Typedef for the iterator type. + */ + typedef typename DoFHandlerType::active_cell_iterator active_cell_iterator; + + /** + * Abbreviation for a pair of + * iterators. + */ + typedef std::pair<active_cell_iterator,active_cell_iterator> iterator_pair; + + /** + * Constructor. Initialize + * the two values by the + * given values. + */ + IteratorRange (const active_cell_iterator &first, + const active_cell_iterator &second); + + /** + * Constructor taking a pair + * of values for + * initialization. + */ + IteratorRange (const iterator_pair &ip); + + /** + * Pair of iterators denoting + * a half-open range. + */ + active_cell_iterator first, second; + }; + + + + + template <typename DoFHandlerType> + inline + IteratorRange<DoFHandlerType>:: + IteratorRange (const active_cell_iterator &first, + const active_cell_iterator &second) + : + first (first), + second (second) + {} + + + + template <typename DoFHandlerType> + inline + IteratorRange<DoFHandlerType>::IteratorRange (const iterator_pair &ip) + : + first (ip.first), + second (ip.second) + {} + + + + namespace AssemblerData + { + template <int dim, + int spacedim> + struct Scratch + { + Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe, + const UpdateFlags update_flags, + const Function<spacedim> *coefficient, + const Function<spacedim> *rhs_function, + const ::dealii::hp::QCollection<dim> &quadrature, + const ::dealii::hp::MappingCollection<dim,spacedim> &mapping) + : + fe_collection (fe), + quadrature_collection (quadrature), + mapping_collection (mapping), + x_fe_values (mapping_collection, + fe_collection, + quadrature_collection, + update_flags), + coefficient_values(quadrature_collection.max_n_quadrature_points()), + coefficient_vector_values (quadrature_collection.max_n_quadrature_points(), + dealii::Vector<double> (fe_collection.n_components())), + rhs_values(quadrature_collection.max_n_quadrature_points()), + rhs_vector_values(quadrature_collection.max_n_quadrature_points(), + dealii::Vector<double> (fe_collection.n_components())), + coefficient (coefficient), + rhs_function (rhs_function), + update_flags (update_flags) + {} + + Scratch (const Scratch &data) + : + fe_collection (data.fe_collection), + quadrature_collection (data.quadrature_collection), + mapping_collection (data.mapping_collection), + x_fe_values (mapping_collection, + fe_collection, + quadrature_collection, + data.update_flags), + coefficient_values (data.coefficient_values), + coefficient_vector_values (data.coefficient_vector_values), + rhs_values (data.rhs_values), + rhs_vector_values (data.rhs_vector_values), + coefficient (data.coefficient), + rhs_function (data.rhs_function), + update_flags (data.update_flags) + {} + + const ::dealii::hp::FECollection<dim,spacedim> &fe_collection; + const ::dealii::hp::QCollection<dim> &quadrature_collection; + const ::dealii::hp::MappingCollection<dim,spacedim> &mapping_collection; + + ::dealii::hp::FEValues<dim,spacedim> x_fe_values; + + std::vector<double> coefficient_values; + std::vector<dealii::Vector<double> > coefficient_vector_values; + std::vector<double> rhs_values; + std::vector<dealii::Vector<double> > rhs_vector_values; + + const Function<spacedim> *coefficient; + const Function<spacedim> *rhs_function; + + const UpdateFlags update_flags; + }; + + + template <typename number> + struct CopyData + { + std::vector<types::global_dof_index> dof_indices; + FullMatrix<number> cell_matrix; + dealii::Vector<number> cell_rhs; + const ConstraintMatrix *constraints; + }; + } + + + template <int dim, + int spacedim, + typename CellIterator, + typename number> + void mass_assembler (const CellIterator &cell, + MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data, + MatrixCreator::internal::AssemblerData::CopyData<number> ©_data) + { + data.x_fe_values.reinit (cell); + const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values (); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + const FiniteElement<dim,spacedim> &fe = fe_values.get_fe(); + const unsigned int n_components = fe.n_components(); + + Assert(data.rhs_function == 0 || + data.rhs_function->n_components==1 || + data.rhs_function->n_components==n_components, + ::dealii::MatrixCreator::ExcComponentMismatch()); + Assert(data.coefficient == 0 || + data.coefficient->n_components==1 || + data.coefficient->n_components==n_components, + ::dealii::MatrixCreator::ExcComponentMismatch()); + + copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell); + copy_data.cell_rhs.reinit (dofs_per_cell); + + copy_data.dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (copy_data.dof_indices); + + const bool use_rhs_function = data.rhs_function != 0; + if (use_rhs_function) + { + if (data.rhs_function->n_components==1) + { + data.rhs_values.resize (n_q_points); + data.rhs_function->value_list (fe_values.get_quadrature_points(), + data.rhs_values); + } + else + { + data.rhs_vector_values.resize (n_q_points, + dealii::Vector<double>(n_components)); + data.rhs_function->vector_value_list (fe_values.get_quadrature_points(), + data.rhs_vector_values); + } + } + + const bool use_coefficient = data.coefficient != 0; + if (use_coefficient) + { + if (data.coefficient->n_components==1) + { + data.coefficient_values.resize (n_q_points); + data.coefficient->value_list (fe_values.get_quadrature_points(), + data.coefficient_values); + } + else + { + data.coefficient_vector_values.resize (n_q_points, + dealii::Vector<double>(n_components)); + data.coefficient->vector_value_list (fe_values.get_quadrature_points(), + data.coefficient_vector_values); + } + } + + + double add_data; + const std::vector<double> &JxW = fe_values.get_JxW_values(); + for (unsigned int i=0; i<dofs_per_cell; ++i) + if (fe.is_primitive ()) + { + const unsigned int component_i = + fe.system_to_component_index(i).first; + const double *phi_i = &fe_values.shape_value(i,0); + add_data = 0; + + // use symmetry in the mass matrix here: + // just need to calculate the diagonal + // and half of the elements above the + // diagonal + for (unsigned int j=i; j<dofs_per_cell; ++j) + if ((n_components==1) || + (fe.system_to_component_index(j).first == + component_i)) + { + const double *phi_j = &fe_values.shape_value(j,0); + add_data = 0; + if (use_coefficient) + { + if (data.coefficient->n_components==1) + for (unsigned int point=0; point<n_q_points; ++point) + add_data += (phi_i[point] * phi_j[point] * JxW[point] * + data.coefficient_values[point]); + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += (phi_i[point] * phi_j[point] * JxW[point] * + data.coefficient_vector_values[point](component_i)); + } + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += phi_i[point] * phi_j[point] * JxW[point]; + + // this is even ok for i==j, since then + // we just write the same value twice. + copy_data.cell_matrix(i,j) = add_data; + copy_data.cell_matrix(j,i) = add_data; + } + + if (use_rhs_function) + { + add_data = 0; + if (data.rhs_function->n_components==1) + for (unsigned int point=0; point<n_q_points; ++point) + add_data += phi_i[point] * JxW[point] * + data.rhs_values[point]; + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += phi_i[point] * JxW[point] * + data.rhs_vector_values[point](component_i); + copy_data.cell_rhs(i) = add_data; + } + } + else + { + // non-primitive vector-valued FE, using + // symmetry again + for (unsigned int j=i; j<dofs_per_cell; ++j) + { + add_data = 0; + for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) + if (fe.get_nonzero_components(i)[comp_i] && + fe.get_nonzero_components(j)[comp_i]) + { + if (use_coefficient) + { + if (data.coefficient->n_components==1) + for (unsigned int point=0; point<n_q_points; ++point) + add_data += (fe_values.shape_value_component(i,point,comp_i) * + fe_values.shape_value_component(j,point,comp_i) * + JxW[point] * + data.coefficient_values[point]); + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += (fe_values.shape_value_component(i,point,comp_i) * + fe_values.shape_value_component(j,point,comp_i) * + JxW[point] * + data.coefficient_vector_values[point](comp_i)); + } + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += fe_values.shape_value_component(i,point,comp_i) * + fe_values.shape_value_component(j,point,comp_i) * JxW[point]; + } + + copy_data.cell_matrix(i,j) = add_data; + copy_data.cell_matrix(j,i) = add_data; + } + + if (use_rhs_function) + { + add_data = 0; + for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) + if (fe.get_nonzero_components(i)[comp_i]) + { + if (data.rhs_function->n_components==1) + for (unsigned int point=0; point<n_q_points; ++point) + add_data += fe_values.shape_value_component(i,point,comp_i) * + JxW[point] * data.rhs_values[point]; + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += fe_values.shape_value_component(i,point,comp_i) * + JxW[point] * data.rhs_vector_values[point](comp_i); + } + copy_data.cell_rhs(i) = add_data; + } + } + } + + + + template <int dim, + int spacedim, + typename CellIterator> + void laplace_assembler (const CellIterator &cell, + MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data, + MatrixCreator::internal::AssemblerData::CopyData<double> ©_data) + { + data.x_fe_values.reinit (cell); + const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values (); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + const FiniteElement<dim,spacedim> &fe = fe_values.get_fe(); + const unsigned int n_components = fe.n_components(); + + Assert(data.rhs_function == 0 || + data.rhs_function->n_components==1 || + data.rhs_function->n_components==n_components, + ::dealii::MatrixCreator::ExcComponentMismatch()); + Assert(data.coefficient == 0 || + data.coefficient->n_components==1 || + data.coefficient->n_components==n_components, + ::dealii::MatrixCreator::ExcComponentMismatch()); + + copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell); + copy_data.cell_rhs.reinit (dofs_per_cell); + copy_data.dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (copy_data.dof_indices); + + + const bool use_rhs_function = data.rhs_function != 0; + if (use_rhs_function) + { + if (data.rhs_function->n_components==1) + { + data.rhs_values.resize (n_q_points); + data.rhs_function->value_list (fe_values.get_quadrature_points(), + data.rhs_values); + } + else + { + data.rhs_vector_values.resize (n_q_points, + dealii::Vector<double>(n_components)); + data.rhs_function->vector_value_list (fe_values.get_quadrature_points(), + data.rhs_vector_values); + } + } + + const bool use_coefficient = data.coefficient != 0; + if (use_coefficient) + { + if (data.coefficient->n_components==1) + { + data.coefficient_values.resize (n_q_points); + data.coefficient->value_list (fe_values.get_quadrature_points(), + data.coefficient_values); + } + else + { + data.coefficient_vector_values.resize (n_q_points, + dealii::Vector<double>(n_components)); + data.coefficient->vector_value_list (fe_values.get_quadrature_points(), + data.coefficient_vector_values); + } + } + + + const std::vector<double> &JxW = fe_values.get_JxW_values(); + double add_data; + for (unsigned int i=0; i<dofs_per_cell; ++i) + if (fe.is_primitive ()) + { + const unsigned int component_i = + fe.system_to_component_index(i).first; + const Tensor<1,spacedim> *grad_phi_i = + &fe_values.shape_grad(i,0); + + // can use symmetry + for (unsigned int j=i; j<dofs_per_cell; ++j) + if ((n_components==1) || + (fe.system_to_component_index(j).first == + component_i)) + { + const Tensor<1,spacedim> *grad_phi_j = + & fe_values.shape_grad(j,0); + add_data = 0; + if (use_coefficient) + { + if (data.coefficient->n_components==1) + for (unsigned int point=0; point<n_q_points; ++point) + add_data += ((grad_phi_i[point]*grad_phi_j[point]) * + JxW[point] * + data.coefficient_values[point]); + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += ((grad_phi_i[point]*grad_phi_j[point]) * + JxW[point] * + data.coefficient_vector_values[point](component_i)); + } + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += (grad_phi_i[point]*grad_phi_j[point]) * + JxW[point]; + + copy_data.cell_matrix(i,j) = add_data; + copy_data.cell_matrix(j,i) = add_data; + } + + if (use_rhs_function) + { + const double *phi_i = &fe_values.shape_value(i,0); + add_data = 0; + if (data.rhs_function->n_components==1) + for (unsigned int point=0; point<n_q_points; ++point) + add_data += phi_i[point] * JxW[point] * + data.rhs_values[point]; + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += phi_i[point] * JxW[point] * + data.rhs_vector_values[point](component_i); + copy_data.cell_rhs(i) = add_data; + } + } + else + { + // non-primitive vector-valued FE + for (unsigned int j=i; j<dofs_per_cell; ++j) + { + add_data = 0; + for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) + if (fe.get_nonzero_components(i)[comp_i] && + fe.get_nonzero_components(j)[comp_i]) + { + if (use_coefficient) + { + if (data.coefficient->n_components==1) + for (unsigned int point=0; point<n_q_points; ++point) + add_data += ((fe_values.shape_grad_component(i,point,comp_i) * + fe_values.shape_grad_component(j,point,comp_i)) * + JxW[point] * + data.coefficient_values[point]); + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += ((fe_values.shape_grad_component(i,point,comp_i) * + fe_values.shape_grad_component(j,point,comp_i)) * + JxW[point] * + data.coefficient_vector_values[point](comp_i)); + } + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += (fe_values.shape_grad_component(i,point,comp_i) * + fe_values.shape_grad_component(j,point,comp_i)) * + JxW[point]; + } + + copy_data.cell_matrix(i,j) = add_data; + copy_data.cell_matrix(j,i) = add_data; + } + + if (use_rhs_function) + { + add_data = 0; + for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) + if (fe.get_nonzero_components(i)[comp_i]) + { + if (data.rhs_function->n_components==1) + for (unsigned int point=0; point<n_q_points; ++point) + add_data += fe_values.shape_value_component(i,point,comp_i) * + JxW[point] * data.rhs_values[point]; + else + for (unsigned int point=0; point<n_q_points; ++point) + add_data += fe_values.shape_value_component(i,point,comp_i) * + JxW[point] * data.rhs_vector_values[point](comp_i); + } + copy_data.cell_rhs(i) = add_data; + } + } + } + + + + template <typename number, + typename MatrixType, + typename VectorType> + void copy_local_to_global (const AssemblerData::CopyData<number> &data, + MatrixType *matrix, + VectorType *right_hand_side) + { + const unsigned int dofs_per_cell = data.dof_indices.size(); + (void)dofs_per_cell; + + Assert (data.cell_matrix.m() == dofs_per_cell, + ExcInternalError()); + Assert (data.cell_matrix.n() == dofs_per_cell, + ExcInternalError()); + Assert ((right_hand_side == 0) + || + (data.cell_rhs.size() == dofs_per_cell), + ExcInternalError()); + + if (right_hand_side != 0) + data.constraints->distribute_local_to_global(data.cell_matrix, + data.cell_rhs, + data.dof_indices, + *matrix, *right_hand_side); + else + data.constraints->distribute_local_to_global(data.cell_matrix, + data.dof_indices, + *matrix); + } + + + + namespace AssemblerBoundary + { + struct Scratch + { + Scratch() {} + }; + + template <typename DoFHandlerType> + struct CopyData + { + CopyData() {}; + + CopyData(CopyData const &data); + + unsigned int dofs_per_cell; + std::vector<types::global_dof_index> dofs; + std::vector<std::vector<bool> > dof_is_on_face; + typename DoFHandlerType::active_cell_iterator cell; + std::vector<FullMatrix<double> > cell_matrix; + std::vector<Vector<double> > cell_vector; + }; + + template <typename DoFHandlerType> + CopyData<DoFHandlerType>::CopyData(CopyData const &data) : + dofs_per_cell(data.dofs_per_cell), + dofs(data.dofs), + dof_is_on_face(data.dof_is_on_face), + cell(data.cell), + cell_matrix(data.cell_matrix), + cell_vector(data.cell_vector) + {} + } + } +} + + +namespace MatrixCreator +{ + + template <int dim, typename number, int spacedim> + void create_mass_matrix (const Mapping<dim,spacedim> &mapping, + const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim> &q, + SparseMatrix<number> &matrix, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + hp::FECollection<dim,spacedim> fe_collection (dof.get_fe()); + hp::QCollection<dim> q_collection (q); + hp::MappingCollection<dim,spacedim> mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> + assembler_data (fe_collection, + update_values | update_JxW_values | + (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), + coefficient, /*rhs_function=*/0, + q_collection, mapping_collection); + + MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.constraints = &constraints; + + WorkStream::run (dof.begin_active(), + static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global<number,SparseMatrix<number>, Vector<number> >, + std_cxx11::_1, &matrix, (Vector<number> *)0), + assembler_data, + copy_data); + } + + + + template <int dim, typename number, int spacedim> + void create_mass_matrix (const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim> &q, + SparseMatrix<number> &matrix, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, + q, matrix, coefficient, constraints); + } + + + + template <int dim, typename number, int spacedim> + void create_mass_matrix (const Mapping<dim,spacedim> &mapping, + const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim> &q, + SparseMatrix<number> &matrix, + const Function<spacedim> &rhs, + Vector<number> &rhs_vector, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + hp::FECollection<dim,spacedim> fe_collection (dof.get_fe()); + hp::QCollection<dim> q_collection (q); + hp::MappingCollection<dim,spacedim> mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> + assembler_data (fe_collection, + update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q_collection, mapping_collection); + MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.constraints = &constraints; + + WorkStream::run (dof.begin_active(), + static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>, + std_cxx11::bind(&MatrixCreator::internal:: + copy_local_to_global<number,SparseMatrix<number>, Vector<number> >, + std_cxx11::_1, &matrix, &rhs_vector), + assembler_data, + copy_data); + } + + + + template <int dim, typename number, int spacedim> + void create_mass_matrix (const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim> &q, + SparseMatrix<number> &matrix, + const Function<spacedim> &rhs, + Vector<number> &rhs_vector, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, + dof, q, matrix, rhs, rhs_vector, coefficient, + constraints); + } + + + + template <int dim, typename number, int spacedim> + void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping, + const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim> &q, + SparseMatrix<number> &matrix, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> + assembler_data (dof.get_fe(), + update_values | update_JxW_values | + (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), + coefficient, /*rhs_function=*/0, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.constraints = &constraints; + + WorkStream::run (dof.begin_active(), + static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global<number,SparseMatrix<number>, Vector<number> >, + std_cxx11::_1, &matrix, (Vector<number> *)0), + assembler_data, + copy_data); + } + + + + template <int dim, typename number, int spacedim> + void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim> &q, + SparseMatrix<number> &matrix, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, + dof, q, matrix, coefficient, constraints); + } + + + + template <int dim, typename number, int spacedim> + void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping, + const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim> &q, + SparseMatrix<number> &matrix, + const Function<spacedim> &rhs, + Vector<number> &rhs_vector, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> + assembler_data (dof.get_fe(), + update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.constraints = &constraints; + + WorkStream::run (dof.begin_active(), + static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global<number,SparseMatrix<number>, Vector<number> >, + std_cxx11::_1, &matrix, &rhs_vector), + assembler_data, + copy_data); + } + + + + template <int dim, typename number, int spacedim> + void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim> &q, + SparseMatrix<number> &matrix, + const Function<spacedim> &rhs, + Vector<number> &rhs_vector, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, + matrix, rhs, rhs_vector, coefficient, constraints); + } + + + + namespace internal + { + template <int dim, int spacedim> + void + static inline + create_boundary_mass_matrix_1 (typename DoFHandler<dim,spacedim>::active_cell_iterator const &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim, + spacedim> > ©_data, + Mapping<dim, spacedim> const &mapping, + FiniteElement<dim,spacedim> const &fe, + Quadrature<dim-1> const &q, + typename FunctionMap<spacedim>::type const &boundary_functions, + Function<spacedim> const *const coefficient, + std::vector<unsigned int> const &component_mapping) + + { + // All assertions for this function + // are in the calling function + // before creating threads. + const unsigned int n_components = fe.n_components(); + const unsigned int n_function_components = boundary_functions.begin()->second->n_components; + const bool fe_is_system = (n_components != 1); + const bool fe_is_primitive = fe.is_primitive(); + + const unsigned int dofs_per_face = fe.dofs_per_face; + + copy_data.cell = cell; + copy_data.dofs_per_cell = fe.dofs_per_cell; + + UpdateFlags update_flags = UpdateFlags (update_values | + update_JxW_values | + update_normal_vectors | + update_quadrature_points); + FEFaceValues<dim,spacedim> fe_values (mapping, fe, q, update_flags); + + // two variables for the coefficient, + // one for the two cases indicated in + // the name + std::vector<double> coefficient_values (fe_values.n_quadrature_points, 1.); + std::vector<Vector<double> > coefficient_vector_values (fe_values.n_quadrature_points, + Vector<double>(n_components)); + const bool coefficient_is_vector = (coefficient != 0 && coefficient->n_components != 1) + ? true : false; + + std::vector<double> rhs_values_scalar (fe_values.n_quadrature_points); + std::vector<Vector<double> > rhs_values_system (fe_values.n_quadrature_points, + Vector<double>(n_function_components)); + + copy_data.dofs.resize(copy_data.dofs_per_cell); + cell->get_dof_indices (copy_data.dofs); + + std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face); + + // Because CopyData objects are reused and that push_back is used, + // dof_is_on_face, cell_matrix, and cell_vector must be cleared before + // they are reused + copy_data.dof_is_on_face.clear(); + copy_data.cell_matrix.clear(); + copy_data.cell_vector.clear(); + + for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(cell->face(face)->boundary_id()) != + boundary_functions.end()) + { + copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell, + copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell)); + fe_values.reinit (cell, face); + + if (fe_is_system) + // FE has several components + { + boundary_functions.find(cell->face(face)->boundary_id()) + ->second->vector_value_list (fe_values.get_quadrature_points(), + rhs_values_system); + + if (coefficient_is_vector) + // If coefficient is + // vector valued, fill + // all components + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_vector_values); + else + { + // If a scalar + // function is + // given, update + // the values, if + // not, use the + // default one set + // in the + // constructor above + if (coefficient != 0) + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + // Copy scalar + // values into vector + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + coefficient_vector_values[point] = coefficient_values[point]; + } + + // Special treatment + // for Hdiv and Hcurl + // elements, where only + // the normal or + // tangential component + // should be projected. + std::vector<std::vector<double> > normal_adjustment(fe_values.n_quadrature_points, + std::vector<double>(n_components, 1.)); + + for (unsigned int comp = 0; comp<n_components; ++comp) + { + const FiniteElement<dim,spacedim> &base = fe.base_element(fe.component_to_base_index(comp).first); + const unsigned int bcomp = fe.component_to_base_index(comp).second; + + if (!base.conforms(FiniteElementData<dim>::H1) && + base.conforms(FiniteElementData<dim>::Hdiv)) + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + normal_adjustment[point][comp] = fe_values.normal_vector(point)[bcomp] + * fe_values.normal_vector(point)[bcomp]; + } + + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + { + const double weight = fe_values.JxW(point); + for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) + if (fe_is_primitive) + { + for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) + { + if (fe.system_to_component_index(j).first + == fe.system_to_component_index(i).first) + { + copy_data.cell_matrix.back()(i,j) + += weight + * fe_values.shape_value(j,point) + * fe_values.shape_value(i,point) + * coefficient_vector_values[point](fe.system_to_component_index(i).first); + } + } + copy_data.cell_vector.back()(i) += fe_values.shape_value(i,point) + * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first]) + * weight; + } + else + { + for (unsigned int comp=0; comp<n_components; ++comp) + { + for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) + copy_data.cell_matrix.back()(i,j) + += fe_values.shape_value_component(j,point,comp) + * fe_values.shape_value_component(i,point,comp) + * normal_adjustment[point][comp] + * weight * coefficient_vector_values[point](comp); + copy_data.cell_vector.back()(i) += fe_values.shape_value_component(i,point,comp) * + rhs_values_system[point](component_mapping[comp]) + * normal_adjustment[point][comp] + * weight; + } + } + } + } + else + // FE is a scalar one + { + boundary_functions.find(cell->face(face)->boundary_id()) + ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); + + if (coefficient != 0) + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + { + const double weight = fe_values.JxW(point); + for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) + { + const double v = fe_values.shape_value(i,point); + for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) + { + const double u = fe_values.shape_value(j,point); + copy_data.cell_matrix.back()(i,j) += (u*v*weight*coefficient_values[point]); + } + copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight; + } + } + } + + + cell->face(face)->get_dof_indices (dofs_on_face_vector); + // for each dof on the cell, have a + // flag whether it is on the face + copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell)); + // check for each of the dofs on this cell + // whether it is on the face + for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) + copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(), + dofs_on_face_vector.end(), + copy_data.dofs[i]) + != + dofs_on_face_vector.end()); + } + } + + template <int dim,int spacedim> + void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim, + spacedim> > const ©_data, + typename FunctionMap<spacedim>::type const &boundary_functions, + std::vector<types::global_dof_index> const &dof_to_boundary_mapping, + SparseMatrix<double> &matrix, + Vector<double> &rhs_vector) + { + // now transfer cell matrix and vector to the whole boundary matrix + // + // in the following: dof[i] holds the global index of the i-th degree of + // freedom on the present cell. If it is also a dof on the boundary, it + // must be a nonzero entry in the dof_to_boundary_mapping and then + // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. + // + // if dof[i] is not on the boundary, it should be zero on the boundary + // therefore on all quadrature points and finally all of its + // entries in the cell matrix and vector should be zero. If not, we + // throw an error (note: because of the evaluation of the shape + // functions only up to machine precision, the term "must be zero" + // really should mean: "should be very small". since this is only an + // assertion and not part of the code, we may choose "very small" + // quite arbitrarily) + // + // the main problem here is that the matrix or vector entry should also + // be zero if the degree of freedom dof[i] is on the boundary, but not + // on the present face, i.e. on another face of the same cell also + // on the boundary. We can therefore not rely on the + // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to + // determine whether dof[i] is a dof on the present face. We do so + // by getting the dofs on the face into @p{dofs_on_face_vector}, + // a vector as always. Usually, searching in a vector is + // inefficient, so we copy the dofs into a set, which enables binary + // searches. + unsigned int pos(0); + for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) + { + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(copy_data.cell->face(face)->boundary_id()) != + boundary_functions.end()) + { + for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) + { + if (copy_data.dof_is_on_face[pos][i] && + dof_to_boundary_mapping[copy_data.dofs[i]] != numbers::invalid_dof_index) + { + for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j) + if (copy_data.dof_is_on_face[pos][j] && + dof_to_boundary_mapping[copy_data.dofs[j]] != numbers::invalid_dof_index) + { + AssertIsFinite(copy_data.cell_matrix[pos](i,j)); + matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]], + dof_to_boundary_mapping[copy_data.dofs[j]], + copy_data.cell_matrix[pos](i,j)); + } + AssertIsFinite(copy_data.cell_vector[pos](i)); + rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += copy_data.cell_vector[pos](i); + } + } + ++pos; + } + } + } + + + template <> + void + inline + create_boundary_mass_matrix_1<1,3> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<1, + 3> > &/*copy_data*/, + Mapping<1,3> const &, + FiniteElement<1,3> const &, + Quadrature<0> const &, + FunctionMap<3>::type const &/*boundary_functions*/, + Function<3> const *const /*coefficient*/, + std::vector<unsigned int> const &/*component_mapping*/) + { + Assert(false,ExcNotImplemented()); + } + } + + + + + template <int dim, int spacedim> + void + create_boundary_mass_matrix (const Mapping<dim, spacedim> &mapping, + const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim-1> &q, + SparseMatrix<double> &matrix, + const typename FunctionMap<spacedim>::type &boundary_functions, + Vector<double> &rhs_vector, + std::vector<types::global_dof_index> &dof_to_boundary_mapping, + const Function<spacedim> *const coefficient, + std::vector<unsigned int> component_mapping) + { + // what would that be in 1d? the + // identity matrix on the boundary + // dofs? + if (dim == 1) + { + Assert (false, ExcNotImplemented()); + return; + } + + const FiniteElement<dim,spacedim> &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components(); + + Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), + ExcInternalError()); + Assert (matrix.n() == matrix.m(), ExcInternalError()); + Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); + Assert (boundary_functions.size() != 0, ExcInternalError()); + Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), + ExcInternalError()); + Assert (coefficient ==0 || + coefficient->n_components==1 || + coefficient->n_components==n_components, ExcComponentMismatch()); + + if (component_mapping.size() == 0) + { + AssertDimension (n_components, boundary_functions.begin()->second->n_components); + for (unsigned int i=0; i<n_components; ++i) + component_mapping.push_back(i); + } + else + AssertDimension (n_components, component_mapping.size()); + + MatrixCreator::internal::AssemblerBoundary::Scratch scratch; + MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast<std_cxx11::function<void (typename DoFHandler<dim,spacedim>::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > &)> > + (std_cxx11::bind(&internal::create_boundary_mass_matrix_1<dim,spacedim>,std_cxx11::_1,std_cxx11::_2, + std_cxx11::_3, + std_cxx11::cref(mapping),std_cxx11::cref(fe),std_cxx11::cref(q), + std_cxx11::cref(boundary_functions),coefficient, + std_cxx11::cref(component_mapping))), + static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary + ::CopyData<DoFHandler<dim,spacedim> > const &)> > (std_cxx11::bind( + &internal::copy_boundary_mass_matrix_1<dim,spacedim>, + std_cxx11::_1, + std_cxx11::cref(boundary_functions), + std_cxx11::cref(dof_to_boundary_mapping), + std_cxx11::ref(matrix), + std_cxx11::ref(rhs_vector))), + scratch, + copy_data); + } + + + + namespace + { + + template <int dim, int spacedim> + void + create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator const + &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary + ::CopyData<hp::DoFHandler<dim,spacedim> > ©_data, + hp::MappingCollection<dim,spacedim> const &mapping, + hp::FECollection<dim,spacedim> const &fe_collection, + hp::QCollection<dim-1> const &q, + const typename FunctionMap<spacedim>::type &boundary_functions, + Function<spacedim> const *const coefficient, + std::vector<unsigned int> const &component_mapping) + { + const unsigned int n_components = fe_collection.n_components(); + const unsigned int n_function_components = boundary_functions.begin()->second->n_components; + const bool fe_is_system = (n_components != 1); + const FiniteElement<dim,spacedim> &fe = cell->get_fe(); + const unsigned int dofs_per_face = fe.dofs_per_face; + + copy_data.cell = cell; + copy_data.dofs_per_cell = fe.dofs_per_cell; + copy_data.dofs.resize(copy_data.dofs_per_cell); + cell->get_dof_indices (copy_data.dofs); + + + UpdateFlags update_flags = UpdateFlags (update_values | + update_JxW_values | + update_quadrature_points); + hp::FEFaceValues<dim,spacedim> x_fe_values (mapping, fe_collection, q, update_flags); + + // two variables for the coefficient, + // one for the two cases indicated in + // the name + std::vector<double> coefficient_values; + std::vector<Vector<double> > coefficient_vector_values; + + std::vector<double> rhs_values_scalar; + std::vector<Vector<double> > rhs_values_system; + + std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face); + + copy_data.dofs.resize(copy_data.dofs_per_cell); + cell->get_dof_indices (copy_data.dofs); + + // Because CopyData objects are reused and that push_back is used, + // dof_is_on_face, cell_matrix, and cell_vector must be cleared before + // they are reused + copy_data.dof_is_on_face.clear(); + copy_data.cell_matrix.clear(); + copy_data.cell_vector.clear(); + + + for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(cell->face(face)->boundary_id()) != + boundary_functions.end()) + { + x_fe_values.reinit (cell, face); + + const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values (); + + copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell, + copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell)); + + if (fe_is_system) + // FE has several components + { + rhs_values_system.resize (fe_values.n_quadrature_points, + Vector<double>(n_function_components)); + boundary_functions.find(cell->face(face)->boundary_id()) + ->second->vector_value_list (fe_values.get_quadrature_points(), + rhs_values_system); + + if (coefficient != 0) + { + if (coefficient->n_components==1) + { + coefficient_values.resize (fe_values.n_quadrature_points); + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + { + const double weight = fe_values.JxW(point); + for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) + { + const double v = fe_values.shape_value(i,point); + for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) + if (fe.system_to_component_index(i).first == + fe.system_to_component_index(j).first) + { + const double u = fe_values.shape_value(j,point); + copy_data.cell_matrix.back()(i,j) + += (u * v * weight * coefficient_values[point]); + } + + copy_data.cell_vector.back()(i) += v * + rhs_values_system[point]( + component_mapping[fe.system_to_component_index(i).first]) * weight; + } + } + } + else + { + coefficient_vector_values.resize (fe_values.n_quadrature_points, + Vector<double>(n_components)); + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_vector_values); + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + { + const double weight = fe_values.JxW(point); + for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) + { + const double v = fe_values.shape_value(i,point); + const unsigned int component_i= + fe.system_to_component_index(i).first; + for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) + if (fe.system_to_component_index(j).first == + component_i) + { + const double u = fe_values.shape_value(j,point); + copy_data.cell_matrix.back()(i,j) += + (u * v * weight * coefficient_vector_values[point](component_i)); + } + copy_data.cell_vector.back()(i) += v * + rhs_values_system[point](component_mapping[component_i]) * weight; + } + } + } + } + else // if (coefficient == 0) + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + { + const double weight = fe_values.JxW(point); + for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) + { + const double v = fe_values.shape_value(i,point); + for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) + if (fe.system_to_component_index(i).first == + fe.system_to_component_index(j).first) + { + const double u = fe_values.shape_value(j,point); + copy_data.cell_matrix.back()(i,j) += (u * v * weight); + } + copy_data.cell_vector.back()(i) += v * + rhs_values_system[point]( + fe.system_to_component_index(i).first) * + weight; + } + } + } + else + // FE is a scalar one + { + rhs_values_scalar.resize (fe_values.n_quadrature_points); + boundary_functions.find(cell->face(face)->boundary_id()) + ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); + + if (coefficient != 0) + { + coefficient_values.resize (fe_values.n_quadrature_points); + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + { + const double weight = fe_values.JxW(point); + for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) + { + const double v = fe_values.shape_value(i,point); + for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) + { + const double u = fe_values.shape_value(j,point); + copy_data.cell_matrix.back()(i,j) += (u * v * weight * + coefficient_values[point]); + } + copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight; + } + } + } + else + for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) + { + const double weight = fe_values.JxW(point); + for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) + { + const double v = fe_values.shape_value(i,point); + for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) + { + const double u = fe_values.shape_value(j,point); + copy_data.cell_matrix.back()(i,j) += (u * v * weight); + } + copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] * weight; + } + } + } + + cell->face(face)->get_dof_indices (dofs_on_face_vector, + cell->active_fe_index()); + // for each dof on the cell, have a + // flag whether it is on the face + copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell)); + // check for each of the dofs on this cell + // whether it is on the face + for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) + copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(), + dofs_on_face_vector.end(), + copy_data.dofs[i]) + != + dofs_on_face_vector.end()); + } + } + + + + template <int dim,int spacedim> + void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary + ::CopyData<hp::DoFHandler<dim,spacedim> > const ©_data, + typename FunctionMap<spacedim>::type const &boundary_functions, + std::vector<types::global_dof_index> const &dof_to_boundary_mapping, + SparseMatrix<double> &matrix, + Vector<double> &rhs_vector) + { + // now transfer cell matrix and vector to the whole boundary matrix + // + // in the following: dof[i] holds the global index of the i-th degree of + // freedom on the present cell. If it is also a dof on the boundary, it + // must be a nonzero entry in the dof_to_boundary_mapping and then + // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. + // + // if dof[i] is not on the boundary, it should be zero on the boundary + // therefore on all quadrature points and finally all of its + // entries in the cell matrix and vector should be zero. If not, we + // throw an error (note: because of the evaluation of the shape + // functions only up to machine precision, the term "must be zero" + // really should mean: "should be very small". since this is only an + // assertion and not part of the code, we may choose "very small" + // quite arbitrarily) + // + // the main problem here is that the matrix or vector entry should also + // be zero if the degree of freedom dof[i] is on the boundary, but not + // on the present face, i.e. on another face of the same cell also + // on the boundary. We can therefore not rely on the + // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to + // determine whether dof[i] is a dof on the present face. We do so + // by getting the dofs on the face into @p{dofs_on_face_vector}, + // a vector as always. Usually, searching in a vector is + // inefficient, so we copy the dofs into a set, which enables binary + // searches. + unsigned int pos(0); + for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) + { + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(copy_data.cell->face(face)->boundary_id()) != + boundary_functions.end()) + { +#ifdef DEBUG + // in debug mode: compute an element in the matrix which is + // guaranteed to belong to a boundary dof. We do this to check that the + // entries in the cell matrix are guaranteed to be zero if the + // respective dof is not on the boundary. Since because of + // round-off, the actual value of the matrix entry may be + // only close to zero, we assert that it is small relative to an element + // which is guaranteed to be nonzero. (absolute smallness does not + // suffice since the size of the domain scales in here) + // + // for this purpose we seek the diagonal of the matrix, where there + // must be an element belonging to the boundary. we take the maximum + // diagonal entry. + types::global_dof_index max_element = static_cast<types::global_dof_index>(0); + for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) && + (*i > max_element)) + max_element = *i; + Assert (max_element == matrix.n()-1, ExcInternalError()); + + double max_diag_entry = 0; + for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) + if (std::fabs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry) + max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i)); +#endif + + for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) + for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j) + { + if (copy_data.dof_is_on_face[pos][i] && copy_data.dof_is_on_face[pos][j]) + matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]], + dof_to_boundary_mapping[copy_data.dofs[j]], + copy_data.cell_matrix[pos](i,j)); + else + { + // assume that all shape functions that are nonzero on the boundary + // are also listed in the @p{dof_to_boundary} mapping. if that + // is not the case, then the boundary mass matrix does not + // make that much sense anyway, as it only contains entries for + // parts of the functions living on the boundary + // + // these, we may compare here for relative smallness of all + // entries in the local matrix which are not taken over to + // the global one + Assert (std::fabs(copy_data.cell_matrix[pos](i,j)) <= 1e-10 * max_diag_entry, + ExcInternalError ()); + } + } + + for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j) + if (copy_data.dof_is_on_face[pos][j]) + rhs_vector(dof_to_boundary_mapping[copy_data.dofs[j]]) += copy_data.cell_vector[pos](j); + else + { + // compare here for relative + // smallness + Assert (std::fabs(copy_data.cell_vector[pos](j)) <= 1e-10 * max_diag_entry, + ExcInternalError()); + } + ++pos; + } + } + } + } + + + + template <int dim, int spacedim> + void create_boundary_mass_matrix (const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim-1> &q, + SparseMatrix<double> &matrix, + const typename FunctionMap<spacedim>::type &rhs, + Vector<double> &rhs_vector, + std::vector<types::global_dof_index> &dof_to_boundary_mapping, + const Function<spacedim> *const a, + std::vector<unsigned int> component_mapping) + { + create_boundary_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q, + matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); + } + + + + template <int dim, int spacedim> + void + create_boundary_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping, + const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim-1> &q, + SparseMatrix<double> &matrix, + const typename FunctionMap<spacedim>::type &boundary_functions, + Vector<double> &rhs_vector, + std::vector<types::global_dof_index> &dof_to_boundary_mapping, + const Function<spacedim> *const coefficient, + std::vector<unsigned int> component_mapping) + { + // what would that be in 1d? the + // identity matrix on the boundary + // dofs? + if (dim == 1) + { + Assert (false, ExcNotImplemented()); + return; + } + + const hp::FECollection<dim,spacedim> &fe_collection = dof.get_fe(); + const unsigned int n_components = fe_collection.n_components(); + + Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), + ExcInternalError()); + Assert (matrix.n() == matrix.m(), ExcInternalError()); + Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); + Assert (boundary_functions.size() != 0, ExcInternalError()); + Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), + ExcInternalError()); + Assert (coefficient ==0 || + coefficient->n_components==1 || + coefficient->n_components==n_components, ExcComponentMismatch()); + + if (component_mapping.size() == 0) + { + AssertDimension (n_components, boundary_functions.begin()->second->n_components); + for (unsigned int i=0; i<n_components; ++i) + component_mapping.push_back(i); + } + else + AssertDimension (n_components, component_mapping.size()); + + MatrixCreator::internal::AssemblerBoundary::Scratch scratch; + MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast<std_cxx11::function<void (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > &)> > + (std_cxx11::bind( &create_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx11::_1,std_cxx11::_2, + std_cxx11::_3, + std_cxx11::cref(mapping),std_cxx11::cref(fe_collection),std_cxx11::cref(q), + std_cxx11::cref(boundary_functions),coefficient, + std_cxx11::cref(component_mapping))), + static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary + ::CopyData<hp::DoFHandler<dim,spacedim> > const &)> > ( + std_cxx11::bind( ©_hp_boundary_mass_matrix_1<dim,spacedim>, + std_cxx11::_1, + std_cxx11::cref(boundary_functions), + std_cxx11::cref(dof_to_boundary_mapping), + std_cxx11::ref(matrix), + std_cxx11::ref(rhs_vector))), + scratch, + copy_data); + } + + + + + template <int dim, int spacedim> + void create_boundary_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim-1> &q, + SparseMatrix<double> &matrix, + const typename FunctionMap<spacedim>::type &rhs, + Vector<double> &rhs_vector, + std::vector<types::global_dof_index> &dof_to_boundary_mapping, + const Function<spacedim> *const a, + std::vector<unsigned int> component_mapping) + { + create_boundary_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, + matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); + } + + + + template <int dim, int spacedim> + void create_laplace_matrix (const Mapping<dim, spacedim> &mapping, + const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim> &q, + SparseMatrix<double> &matrix, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + hp::FECollection<dim,spacedim> fe_collection (dof.get_fe()); + hp::QCollection<dim> q_collection (q); + hp::MappingCollection<dim,spacedim> mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> + assembler_data (fe_collection, + update_gradients | update_JxW_values | + (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), + coefficient, /*rhs_function=*/0, + q_collection, mapping_collection); + MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.constraints = &constraints; + + WorkStream::run (dof.begin_active(), + static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global<double,SparseMatrix<double>, Vector<double> >, + std_cxx11::_1, + &matrix, + (Vector<double> *)NULL), + assembler_data, + copy_data); + } + + + + template <int dim, int spacedim> + void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim> &q, + SparseMatrix<double> &matrix, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, + dof, q, matrix, coefficient, constraints); + } + + + + template <int dim, int spacedim> + void create_laplace_matrix (const Mapping<dim, spacedim> &mapping, + const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim> &q, + SparseMatrix<double> &matrix, + const Function<spacedim> &rhs, + Vector<double> &rhs_vector, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + hp::FECollection<dim,spacedim> fe_collection (dof.get_fe()); + hp::QCollection<dim> q_collection (q); + hp::MappingCollection<dim,spacedim> mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> + assembler_data (fe_collection, + update_gradients | update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q_collection, mapping_collection); + MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.constraints = &constraints; + + WorkStream::run (dof.begin_active(), + static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global<double,SparseMatrix<double>, Vector<double> >, + std_cxx11::_1, + &matrix, + &rhs_vector), + assembler_data, + copy_data); + } + + + + template <int dim, int spacedim> + void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof, + const Quadrature<dim> &q, + SparseMatrix<double> &matrix, + const Function<spacedim> &rhs, + Vector<double> &rhs_vector, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q, + matrix, rhs, rhs_vector, coefficient, constraints); + } + + + + template <int dim, int spacedim> + void create_laplace_matrix (const hp::MappingCollection<dim,spacedim> &mapping, + const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim> &q, + SparseMatrix<double> &matrix, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> + assembler_data (dof.get_fe(), + update_gradients | update_JxW_values | + (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), + coefficient, /*rhs_function=*/0, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.constraints = &constraints; + + WorkStream::run (dof.begin_active(), + static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global<double,SparseMatrix<double>, Vector<double> >, + std_cxx11::_1, + &matrix, + (Vector<double> *)0), + assembler_data, + copy_data); + } + + + + template <int dim, int spacedim> + void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim> &q, + SparseMatrix<double> &matrix, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, + matrix, coefficient, constraints); + } + + + + template <int dim, int spacedim> + void create_laplace_matrix (const hp::MappingCollection<dim,spacedim> &mapping, + const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim> &q, + SparseMatrix<double> &matrix, + const Function<spacedim> &rhs, + Vector<double> &rhs_vector, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> + assembler_data (dof.get_fe(), + update_gradients | update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.constraints = &constraints; + + WorkStream::run (dof.begin_active(), + static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global<double,SparseMatrix<double>, Vector<double> >, + std_cxx11::_1, + &matrix, + &rhs_vector), + assembler_data, + copy_data); + } + + + + template <int dim, int spacedim> + void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof, + const hp::QCollection<dim> &q, + SparseMatrix<double> &matrix, + const Function<spacedim> &rhs, + Vector<double> &rhs_vector, + const Function<spacedim> *const coefficient, + const ConstraintMatrix &constraints) + { + create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, + matrix, rhs, rhs_vector, coefficient, constraints); + } + +} // namespace MatrixCreator + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/numerics/CMakeLists.txt b/source/numerics/CMakeLists.txt index 086368040c..c8adfa339b 100644 --- a/source/numerics/CMakeLists.txt +++ b/source/numerics/CMakeLists.txt @@ -29,6 +29,10 @@ SET(_src error_estimator_inst2.cc fe_field_function.cc histogram.cc + matrix_creator.cc + matrix_creator_inst2.cc + matrix_creator_inst3.cc + matrix_tools_once.cc matrix_tools.cc point_value_history.cc solution_transfer.cc @@ -59,6 +63,7 @@ SET(_inst error_estimator_1d.inst.in error_estimator.inst.in fe_field_function.inst.in + matrix_creator.inst.in matrix_tools.inst.in point_value_history.inst.in solution_transfer.inst.in diff --git a/source/numerics/matrix_creator.cc b/source/numerics/matrix_creator.cc new file mode 100644 index 0000000000..1bb6205054 --- /dev/null +++ b/source/numerics/matrix_creator.cc @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2015 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. +// +// --------------------------------------------------------------------- + +#include <deal.II/numerics/matrix_creator.templates.h> + + +DEAL_II_NAMESPACE_OPEN + +// explicit instantiations +#define SPLIT_INSTANTIATIONS_COUNT 3 +#ifndef SPLIT_INSTANTIATIONS_INDEX +#define SPLIT_INSTANTIATIONS_INDEX 0 +#endif +#include "matrix_creator.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/numerics/matrix_creator.inst.in b/source/numerics/matrix_creator.inst.in new file mode 100644 index 0000000000..a84cafccb1 --- /dev/null +++ b/source/numerics/matrix_creator.inst.in @@ -0,0 +1,315 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2015 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. +// +// --------------------------------------------------------------------- + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + +// non-hp version of create_mass_matrix + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> + (const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping, + const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> + (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> + (const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping, + const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_space_dimension> &rhs, + Vector<double> &rhs_vector, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> + (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_space_dimension> &rhs, + Vector<double> &rhs_vector, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + + + template + void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension> + (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension-1> &q, + SparseMatrix<double> &matrix, + const FunctionMap<deal_II_space_dimension>::type &rhs, + Vector<double> &rhs_vector, + std::vector<types::global_dof_index> &dof_to_boundary_mapping, + const Function<deal_II_space_dimension> * const a, + std::vector<unsigned int>); + + template + void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension> + (const Mapping<deal_II_dimension,deal_II_space_dimension> &, + const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension-1> &q, + SparseMatrix<double> &matrix, + const FunctionMap<deal_II_space_dimension>::type &rhs, + Vector<double> &rhs_vector, + std::vector<types::global_dof_index> &dof_to_boundary_mapping, + const Function<deal_II_space_dimension> * const a, + std::vector<unsigned int>); + + template + void + MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension> + (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>&, + const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&, + const hp::QCollection<deal_II_dimension-1>&, + SparseMatrix<double>&, + const FunctionMap<deal_II_space_dimension>::type&, + Vector<double>&, + std::vector<types::global_dof_index>&, + const Function<deal_II_space_dimension> * const, + std::vector<unsigned int>); + + template + void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension> + (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&, + const hp::QCollection<deal_II_dimension-1>&, + SparseMatrix<double>&, + const FunctionMap<deal_II_space_dimension>::type&, + Vector<double>&, + std::vector<types::global_dof_index>&, + const Function<deal_II_space_dimension> * const, + std::vector<unsigned int>); + +// same for float + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> + (const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping, + const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<float> &matrix, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> + (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<float> &matrix, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> + (const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping, + const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<float> &matrix, + const Function<deal_II_space_dimension> &rhs, + Vector<float> &rhs_vector, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> + (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<float> &matrix, + const Function<deal_II_space_dimension> &rhs, + Vector<float> &rhs_vector, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + +#endif + } + + +//TODO[SP]: replace <deal_II_dimension> by <deal_II_dimension, deal_II_space_dimension> +// where applicable and move to codimension cases above also when applicable +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +// hp versions of functions +#if deal_II_dimension <= deal_II_space_dimension + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> + (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping, + const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> + (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping, + const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_space_dimension> &rhs, + Vector<double> &rhs_vector, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + + +#endif + +#if deal_II_dimension == deal_II_space_dimension + + + template + void MatrixCreator::create_mass_matrix<deal_II_dimension> + (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix<deal_II_dimension> + (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_space_dimension> &rhs, + Vector<double> &rhs_vector, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + +// non-hp versions of create_laplace_matrix + template + void MatrixCreator::create_laplace_matrix<deal_II_dimension> + (const DoFHandler<deal_II_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix<deal_II_dimension> + (const Mapping<deal_II_dimension> &mapping, + const DoFHandler<deal_II_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix<deal_II_dimension> + (const Mapping<deal_II_dimension> &mapping, + const DoFHandler<deal_II_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_dimension> &rhs, + Vector<double> &rhs_vector, + const Function<deal_II_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix<deal_II_dimension> + (const DoFHandler<deal_II_dimension> &dof, + const Quadrature<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_dimension> &rhs, + Vector<double> &rhs_vector, + const Function<deal_II_dimension> * const coefficient, + const ConstraintMatrix &constraints); + +// hp versions of create_laplace_matrix + template + void MatrixCreator::create_laplace_matrix<deal_II_dimension> + (const hp::DoFHandler<deal_II_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix<deal_II_dimension> + (const hp::MappingCollection<deal_II_dimension> &mapping, + const hp::DoFHandler<deal_II_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix<deal_II_dimension> + (const hp::MappingCollection<deal_II_dimension> &mapping, + const hp::DoFHandler<deal_II_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_dimension> &rhs, + Vector<double> &rhs_vector, + const Function<deal_II_dimension> * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix<deal_II_dimension> + (const hp::DoFHandler<deal_II_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<double> &matrix, + const Function<deal_II_dimension> &rhs, + Vector<double> &rhs_vector, + const Function<deal_II_dimension> * const coefficient, + const ConstraintMatrix &constraints); + +#endif + +// same for float +#if deal_II_dimension <= deal_II_space_dimension + + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> + (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping, + const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<float> &matrix, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> + (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping, + const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<float> &matrix, + const Function<deal_II_space_dimension> &rhs, + Vector<float> &rhs_vector, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + + +#endif + +#if deal_II_dimension == deal_II_space_dimension + + + template + void MatrixCreator::create_mass_matrix<deal_II_dimension> + (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<float> &matrix, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix<deal_II_dimension> + (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, + const hp::QCollection<deal_II_dimension> &q, + SparseMatrix<float> &matrix, + const Function<deal_II_space_dimension> &rhs, + Vector<float> &rhs_vector, + const Function<deal_II_space_dimension> * const coefficient, + const ConstraintMatrix &constraints); + +#endif diff --git a/source/numerics/matrix_creator_inst2.cc b/source/numerics/matrix_creator_inst2.cc new file mode 100644 index 0000000000..73c34c503b --- /dev/null +++ b/source/numerics/matrix_creator_inst2.cc @@ -0,0 +1,20 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + +// This file compiles a part of the instantiations from matrix_creator.cc +// to reduce the compilation unit (and memory consumption) + +#define SPLIT_INSTANTIATIONS_INDEX 1 +#include "matrix_creator.cc" diff --git a/source/numerics/matrix_creator_inst3.cc b/source/numerics/matrix_creator_inst3.cc new file mode 100644 index 0000000000..2cc0df2370 --- /dev/null +++ b/source/numerics/matrix_creator_inst3.cc @@ -0,0 +1,20 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + +// This file compiles a part of the instantiations from matrix_creator.cc +// to reduce the compilation unit (and memory consumption) + +#define SPLIT_INSTANTIATIONS_INDEX 2 +#include "matrix_creator.cc" diff --git a/source/numerics/matrix_tools.cc b/source/numerics/matrix_tools.cc index ddd5358422..28b46756c1 100644 --- a/source/numerics/matrix_tools.cc +++ b/source/numerics/matrix_tools.cc @@ -59,1833 +59,6 @@ DEAL_II_NAMESPACE_OPEN - -namespace MatrixCreator -{ - namespace internal - { - /** - * Convenience abbreviation for - * pairs of DoF handler cell - * iterators. This type works - * just like a - * <tt>std::pair<iterator,iterator></tt> - * but is templatized on the - * dof handler that should be used. - */ - template <typename DoFHandlerType> - struct IteratorRange - { - /** - * Typedef for the iterator type. - */ - typedef typename DoFHandlerType::active_cell_iterator active_cell_iterator; - - /** - * Abbreviation for a pair of - * iterators. - */ - typedef std::pair<active_cell_iterator,active_cell_iterator> iterator_pair; - - /** - * Constructor. Initialize - * the two values by the - * given values. - */ - IteratorRange (const active_cell_iterator &first, - const active_cell_iterator &second); - - /** - * Constructor taking a pair - * of values for - * initialization. - */ - IteratorRange (const iterator_pair &ip); - - /** - * Pair of iterators denoting - * a half-open range. - */ - active_cell_iterator first, second; - }; - - - - - template <typename DoFHandlerType> - inline - IteratorRange<DoFHandlerType>:: - IteratorRange (const active_cell_iterator &first, - const active_cell_iterator &second) - : - first (first), - second (second) - {} - - - - template <typename DoFHandlerType> - inline - IteratorRange<DoFHandlerType>::IteratorRange (const iterator_pair &ip) - : - first (ip.first), - second (ip.second) - {} - - - - namespace AssemblerData - { - template <int dim, - int spacedim> - struct Scratch - { - Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe, - const UpdateFlags update_flags, - const Function<spacedim> *coefficient, - const Function<spacedim> *rhs_function, - const ::dealii::hp::QCollection<dim> &quadrature, - const ::dealii::hp::MappingCollection<dim,spacedim> &mapping) - : - fe_collection (fe), - quadrature_collection (quadrature), - mapping_collection (mapping), - x_fe_values (mapping_collection, - fe_collection, - quadrature_collection, - update_flags), - coefficient_values(quadrature_collection.max_n_quadrature_points()), - coefficient_vector_values (quadrature_collection.max_n_quadrature_points(), - dealii::Vector<double> (fe_collection.n_components())), - rhs_values(quadrature_collection.max_n_quadrature_points()), - rhs_vector_values(quadrature_collection.max_n_quadrature_points(), - dealii::Vector<double> (fe_collection.n_components())), - coefficient (coefficient), - rhs_function (rhs_function), - update_flags (update_flags) - {} - - Scratch (const Scratch &data) - : - fe_collection (data.fe_collection), - quadrature_collection (data.quadrature_collection), - mapping_collection (data.mapping_collection), - x_fe_values (mapping_collection, - fe_collection, - quadrature_collection, - data.update_flags), - coefficient_values (data.coefficient_values), - coefficient_vector_values (data.coefficient_vector_values), - rhs_values (data.rhs_values), - rhs_vector_values (data.rhs_vector_values), - coefficient (data.coefficient), - rhs_function (data.rhs_function), - update_flags (data.update_flags) - {} - - const ::dealii::hp::FECollection<dim,spacedim> &fe_collection; - const ::dealii::hp::QCollection<dim> &quadrature_collection; - const ::dealii::hp::MappingCollection<dim,spacedim> &mapping_collection; - - ::dealii::hp::FEValues<dim,spacedim> x_fe_values; - - std::vector<double> coefficient_values; - std::vector<dealii::Vector<double> > coefficient_vector_values; - std::vector<double> rhs_values; - std::vector<dealii::Vector<double> > rhs_vector_values; - - const Function<spacedim> *coefficient; - const Function<spacedim> *rhs_function; - - const UpdateFlags update_flags; - }; - - - template <typename number> - struct CopyData - { - std::vector<types::global_dof_index> dof_indices; - FullMatrix<number> cell_matrix; - dealii::Vector<number> cell_rhs; - const ConstraintMatrix *constraints; - }; - } - - - template <int dim, - int spacedim, - typename CellIterator, - typename number> - void mass_assembler (const CellIterator &cell, - MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data, - MatrixCreator::internal::AssemblerData::CopyData<number> ©_data) - { - data.x_fe_values.reinit (cell); - const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values (); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - const FiniteElement<dim,spacedim> &fe = fe_values.get_fe(); - const unsigned int n_components = fe.n_components(); - - Assert(data.rhs_function == 0 || - data.rhs_function->n_components==1 || - data.rhs_function->n_components==n_components, - ::dealii::MatrixCreator::ExcComponentMismatch()); - Assert(data.coefficient == 0 || - data.coefficient->n_components==1 || - data.coefficient->n_components==n_components, - ::dealii::MatrixCreator::ExcComponentMismatch()); - - copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell); - copy_data.cell_rhs.reinit (dofs_per_cell); - - copy_data.dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (copy_data.dof_indices); - - const bool use_rhs_function = data.rhs_function != 0; - if (use_rhs_function) - { - if (data.rhs_function->n_components==1) - { - data.rhs_values.resize (n_q_points); - data.rhs_function->value_list (fe_values.get_quadrature_points(), - data.rhs_values); - } - else - { - data.rhs_vector_values.resize (n_q_points, - dealii::Vector<double>(n_components)); - data.rhs_function->vector_value_list (fe_values.get_quadrature_points(), - data.rhs_vector_values); - } - } - - const bool use_coefficient = data.coefficient != 0; - if (use_coefficient) - { - if (data.coefficient->n_components==1) - { - data.coefficient_values.resize (n_q_points); - data.coefficient->value_list (fe_values.get_quadrature_points(), - data.coefficient_values); - } - else - { - data.coefficient_vector_values.resize (n_q_points, - dealii::Vector<double>(n_components)); - data.coefficient->vector_value_list (fe_values.get_quadrature_points(), - data.coefficient_vector_values); - } - } - - - double add_data; - const std::vector<double> &JxW = fe_values.get_JxW_values(); - for (unsigned int i=0; i<dofs_per_cell; ++i) - if (fe.is_primitive ()) - { - const unsigned int component_i = - fe.system_to_component_index(i).first; - const double *phi_i = &fe_values.shape_value(i,0); - add_data = 0; - - // use symmetry in the mass matrix here: - // just need to calculate the diagonal - // and half of the elements above the - // diagonal - for (unsigned int j=i; j<dofs_per_cell; ++j) - if ((n_components==1) || - (fe.system_to_component_index(j).first == - component_i)) - { - const double *phi_j = &fe_values.shape_value(j,0); - add_data = 0; - if (use_coefficient) - { - if (data.coefficient->n_components==1) - for (unsigned int point=0; point<n_q_points; ++point) - add_data += (phi_i[point] * phi_j[point] * JxW[point] * - data.coefficient_values[point]); - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += (phi_i[point] * phi_j[point] * JxW[point] * - data.coefficient_vector_values[point](component_i)); - } - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += phi_i[point] * phi_j[point] * JxW[point]; - - // this is even ok for i==j, since then - // we just write the same value twice. - copy_data.cell_matrix(i,j) = add_data; - copy_data.cell_matrix(j,i) = add_data; - } - - if (use_rhs_function) - { - add_data = 0; - if (data.rhs_function->n_components==1) - for (unsigned int point=0; point<n_q_points; ++point) - add_data += phi_i[point] * JxW[point] * - data.rhs_values[point]; - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += phi_i[point] * JxW[point] * - data.rhs_vector_values[point](component_i); - copy_data.cell_rhs(i) = add_data; - } - } - else - { - // non-primitive vector-valued FE, using - // symmetry again - for (unsigned int j=i; j<dofs_per_cell; ++j) - { - add_data = 0; - for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) - if (fe.get_nonzero_components(i)[comp_i] && - fe.get_nonzero_components(j)[comp_i]) - { - if (use_coefficient) - { - if (data.coefficient->n_components==1) - for (unsigned int point=0; point<n_q_points; ++point) - add_data += (fe_values.shape_value_component(i,point,comp_i) * - fe_values.shape_value_component(j,point,comp_i) * - JxW[point] * - data.coefficient_values[point]); - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += (fe_values.shape_value_component(i,point,comp_i) * - fe_values.shape_value_component(j,point,comp_i) * - JxW[point] * - data.coefficient_vector_values[point](comp_i)); - } - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += fe_values.shape_value_component(i,point,comp_i) * - fe_values.shape_value_component(j,point,comp_i) * JxW[point]; - } - - copy_data.cell_matrix(i,j) = add_data; - copy_data.cell_matrix(j,i) = add_data; - } - - if (use_rhs_function) - { - add_data = 0; - for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) - if (fe.get_nonzero_components(i)[comp_i]) - { - if (data.rhs_function->n_components==1) - for (unsigned int point=0; point<n_q_points; ++point) - add_data += fe_values.shape_value_component(i,point,comp_i) * - JxW[point] * data.rhs_values[point]; - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += fe_values.shape_value_component(i,point,comp_i) * - JxW[point] * data.rhs_vector_values[point](comp_i); - } - copy_data.cell_rhs(i) = add_data; - } - } - } - - - - template <int dim, - int spacedim, - typename CellIterator> - void laplace_assembler (const CellIterator &cell, - MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data, - MatrixCreator::internal::AssemblerData::CopyData<double> ©_data) - { - data.x_fe_values.reinit (cell); - const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values (); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - const FiniteElement<dim,spacedim> &fe = fe_values.get_fe(); - const unsigned int n_components = fe.n_components(); - - Assert(data.rhs_function == 0 || - data.rhs_function->n_components==1 || - data.rhs_function->n_components==n_components, - ::dealii::MatrixCreator::ExcComponentMismatch()); - Assert(data.coefficient == 0 || - data.coefficient->n_components==1 || - data.coefficient->n_components==n_components, - ::dealii::MatrixCreator::ExcComponentMismatch()); - - copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell); - copy_data.cell_rhs.reinit (dofs_per_cell); - copy_data.dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (copy_data.dof_indices); - - - const bool use_rhs_function = data.rhs_function != 0; - if (use_rhs_function) - { - if (data.rhs_function->n_components==1) - { - data.rhs_values.resize (n_q_points); - data.rhs_function->value_list (fe_values.get_quadrature_points(), - data.rhs_values); - } - else - { - data.rhs_vector_values.resize (n_q_points, - dealii::Vector<double>(n_components)); - data.rhs_function->vector_value_list (fe_values.get_quadrature_points(), - data.rhs_vector_values); - } - } - - const bool use_coefficient = data.coefficient != 0; - if (use_coefficient) - { - if (data.coefficient->n_components==1) - { - data.coefficient_values.resize (n_q_points); - data.coefficient->value_list (fe_values.get_quadrature_points(), - data.coefficient_values); - } - else - { - data.coefficient_vector_values.resize (n_q_points, - dealii::Vector<double>(n_components)); - data.coefficient->vector_value_list (fe_values.get_quadrature_points(), - data.coefficient_vector_values); - } - } - - - const std::vector<double> &JxW = fe_values.get_JxW_values(); - double add_data; - for (unsigned int i=0; i<dofs_per_cell; ++i) - if (fe.is_primitive ()) - { - const unsigned int component_i = - fe.system_to_component_index(i).first; - const Tensor<1,spacedim> *grad_phi_i = - &fe_values.shape_grad(i,0); - - // can use symmetry - for (unsigned int j=i; j<dofs_per_cell; ++j) - if ((n_components==1) || - (fe.system_to_component_index(j).first == - component_i)) - { - const Tensor<1,spacedim> *grad_phi_j = - & fe_values.shape_grad(j,0); - add_data = 0; - if (use_coefficient) - { - if (data.coefficient->n_components==1) - for (unsigned int point=0; point<n_q_points; ++point) - add_data += ((grad_phi_i[point]*grad_phi_j[point]) * - JxW[point] * - data.coefficient_values[point]); - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += ((grad_phi_i[point]*grad_phi_j[point]) * - JxW[point] * - data.coefficient_vector_values[point](component_i)); - } - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += (grad_phi_i[point]*grad_phi_j[point]) * - JxW[point]; - - copy_data.cell_matrix(i,j) = add_data; - copy_data.cell_matrix(j,i) = add_data; - } - - if (use_rhs_function) - { - const double *phi_i = &fe_values.shape_value(i,0); - add_data = 0; - if (data.rhs_function->n_components==1) - for (unsigned int point=0; point<n_q_points; ++point) - add_data += phi_i[point] * JxW[point] * - data.rhs_values[point]; - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += phi_i[point] * JxW[point] * - data.rhs_vector_values[point](component_i); - copy_data.cell_rhs(i) = add_data; - } - } - else - { - // non-primitive vector-valued FE - for (unsigned int j=i; j<dofs_per_cell; ++j) - { - add_data = 0; - for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) - if (fe.get_nonzero_components(i)[comp_i] && - fe.get_nonzero_components(j)[comp_i]) - { - if (use_coefficient) - { - if (data.coefficient->n_components==1) - for (unsigned int point=0; point<n_q_points; ++point) - add_data += ((fe_values.shape_grad_component(i,point,comp_i) * - fe_values.shape_grad_component(j,point,comp_i)) * - JxW[point] * - data.coefficient_values[point]); - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += ((fe_values.shape_grad_component(i,point,comp_i) * - fe_values.shape_grad_component(j,point,comp_i)) * - JxW[point] * - data.coefficient_vector_values[point](comp_i)); - } - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += (fe_values.shape_grad_component(i,point,comp_i) * - fe_values.shape_grad_component(j,point,comp_i)) * - JxW[point]; - } - - copy_data.cell_matrix(i,j) = add_data; - copy_data.cell_matrix(j,i) = add_data; - } - - if (use_rhs_function) - { - add_data = 0; - for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) - if (fe.get_nonzero_components(i)[comp_i]) - { - if (data.rhs_function->n_components==1) - for (unsigned int point=0; point<n_q_points; ++point) - add_data += fe_values.shape_value_component(i,point,comp_i) * - JxW[point] * data.rhs_values[point]; - else - for (unsigned int point=0; point<n_q_points; ++point) - add_data += fe_values.shape_value_component(i,point,comp_i) * - JxW[point] * data.rhs_vector_values[point](comp_i); - } - copy_data.cell_rhs(i) = add_data; - } - } - } - - - - template <typename number, - typename MatrixType, - typename VectorType> - void copy_local_to_global (const AssemblerData::CopyData<number> &data, - MatrixType *matrix, - VectorType *right_hand_side) - { - const unsigned int dofs_per_cell = data.dof_indices.size(); - (void)dofs_per_cell; - - Assert (data.cell_matrix.m() == dofs_per_cell, - ExcInternalError()); - Assert (data.cell_matrix.n() == dofs_per_cell, - ExcInternalError()); - Assert ((right_hand_side == 0) - || - (data.cell_rhs.size() == dofs_per_cell), - ExcInternalError()); - - if (right_hand_side != 0) - data.constraints->distribute_local_to_global(data.cell_matrix, - data.cell_rhs, - data.dof_indices, - *matrix, *right_hand_side); - else - data.constraints->distribute_local_to_global(data.cell_matrix, - data.dof_indices, - *matrix); - } - - - - namespace AssemblerBoundary - { - struct Scratch - { - Scratch() {} - }; - - template <typename DoFHandlerType> - struct CopyData - { - CopyData() {}; - - CopyData(CopyData const &data); - - unsigned int dofs_per_cell; - std::vector<types::global_dof_index> dofs; - std::vector<std::vector<bool> > dof_is_on_face; - typename DoFHandlerType::active_cell_iterator cell; - std::vector<FullMatrix<double> > cell_matrix; - std::vector<Vector<double> > cell_vector; - }; - - template <typename DoFHandlerType> - CopyData<DoFHandlerType>::CopyData(CopyData const &data) : - dofs_per_cell(data.dofs_per_cell), - dofs(data.dofs), - dof_is_on_face(data.dof_is_on_face), - cell(data.cell), - cell_matrix(data.cell_matrix), - cell_vector(data.cell_vector) - {} - } - } -} - - -namespace MatrixCreator -{ - - template <int dim, typename number, int spacedim> - void create_mass_matrix (const Mapping<dim,spacedim> &mapping, - const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim> &q, - SparseMatrix<number> &matrix, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - hp::FECollection<dim,spacedim> fe_collection (dof.get_fe()); - hp::QCollection<dim> q_collection (q); - hp::MappingCollection<dim,spacedim> mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> - assembler_data (fe_collection, - update_values | update_JxW_values | - (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), - coefficient, /*rhs_function=*/0, - q_collection, mapping_collection); - - MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; - - WorkStream::run (dof.begin_active(), - static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global<number,SparseMatrix<number>, Vector<number> >, - std_cxx11::_1, &matrix, (Vector<number> *)0), - assembler_data, - copy_data); - } - - - - template <int dim, typename number, int spacedim> - void create_mass_matrix (const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim> &q, - SparseMatrix<number> &matrix, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, - q, matrix, coefficient, constraints); - } - - - - template <int dim, typename number, int spacedim> - void create_mass_matrix (const Mapping<dim,spacedim> &mapping, - const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim> &q, - SparseMatrix<number> &matrix, - const Function<spacedim> &rhs, - Vector<number> &rhs_vector, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - hp::FECollection<dim,spacedim> fe_collection (dof.get_fe()); - hp::QCollection<dim> q_collection (q); - hp::MappingCollection<dim,spacedim> mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> - assembler_data (fe_collection, - update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q_collection, mapping_collection); - MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; - - WorkStream::run (dof.begin_active(), - static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>, - std_cxx11::bind(&MatrixCreator::internal:: - copy_local_to_global<number,SparseMatrix<number>, Vector<number> >, - std_cxx11::_1, &matrix, &rhs_vector), - assembler_data, - copy_data); - } - - - - template <int dim, typename number, int spacedim> - void create_mass_matrix (const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim> &q, - SparseMatrix<number> &matrix, - const Function<spacedim> &rhs, - Vector<number> &rhs_vector, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, - dof, q, matrix, rhs, rhs_vector, coefficient, - constraints); - } - - - - template <int dim, typename number, int spacedim> - void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping, - const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim> &q, - SparseMatrix<number> &matrix, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> - assembler_data (dof.get_fe(), - update_values | update_JxW_values | - (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), - coefficient, /*rhs_function=*/0, - q, mapping); - MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; - - WorkStream::run (dof.begin_active(), - static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global<number,SparseMatrix<number>, Vector<number> >, - std_cxx11::_1, &matrix, (Vector<number> *)0), - assembler_data, - copy_data); - } - - - - template <int dim, typename number, int spacedim> - void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim> &q, - SparseMatrix<number> &matrix, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, - dof, q, matrix, coefficient, constraints); - } - - - - template <int dim, typename number, int spacedim> - void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping, - const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim> &q, - SparseMatrix<number> &matrix, - const Function<spacedim> &rhs, - Vector<number> &rhs_vector, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> - assembler_data (dof.get_fe(), - update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q, mapping); - MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; - - WorkStream::run (dof.begin_active(), - static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global<number,SparseMatrix<number>, Vector<number> >, - std_cxx11::_1, &matrix, &rhs_vector), - assembler_data, - copy_data); - } - - - - template <int dim, typename number, int spacedim> - void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim> &q, - SparseMatrix<number> &matrix, - const Function<spacedim> &rhs, - Vector<number> &rhs_vector, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, - matrix, rhs, rhs_vector, coefficient, constraints); - } - - - - namespace - { - template <int dim, int spacedim> - void - create_boundary_mass_matrix_1 (typename DoFHandler<dim,spacedim>::active_cell_iterator const &cell, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim, - spacedim> > ©_data, - Mapping<dim, spacedim> const &mapping, - FiniteElement<dim,spacedim> const &fe, - Quadrature<dim-1> const &q, - typename FunctionMap<spacedim>::type const &boundary_functions, - Function<spacedim> const *const coefficient, - std::vector<unsigned int> const &component_mapping) - - { - // All assertions for this function - // are in the calling function - // before creating threads. - const unsigned int n_components = fe.n_components(); - const unsigned int n_function_components = boundary_functions.begin()->second->n_components; - const bool fe_is_system = (n_components != 1); - const bool fe_is_primitive = fe.is_primitive(); - - const unsigned int dofs_per_face = fe.dofs_per_face; - - copy_data.cell = cell; - copy_data.dofs_per_cell = fe.dofs_per_cell; - - UpdateFlags update_flags = UpdateFlags (update_values | - update_JxW_values | - update_normal_vectors | - update_quadrature_points); - FEFaceValues<dim,spacedim> fe_values (mapping, fe, q, update_flags); - - // two variables for the coefficient, - // one for the two cases indicated in - // the name - std::vector<double> coefficient_values (fe_values.n_quadrature_points, 1.); - std::vector<Vector<double> > coefficient_vector_values (fe_values.n_quadrature_points, - Vector<double>(n_components)); - const bool coefficient_is_vector = (coefficient != 0 && coefficient->n_components != 1) - ? true : false; - - std::vector<double> rhs_values_scalar (fe_values.n_quadrature_points); - std::vector<Vector<double> > rhs_values_system (fe_values.n_quadrature_points, - Vector<double>(n_function_components)); - - copy_data.dofs.resize(copy_data.dofs_per_cell); - cell->get_dof_indices (copy_data.dofs); - - std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face); - - // Because CopyData objects are reused and that push_back is used, - // dof_is_on_face, cell_matrix, and cell_vector must be cleared before - // they are reused - copy_data.dof_is_on_face.clear(); - copy_data.cell_matrix.clear(); - copy_data.cell_vector.clear(); - - for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) - // check if this face is on that part of - // the boundary we are interested in - if (boundary_functions.find(cell->face(face)->boundary_id()) != - boundary_functions.end()) - { - copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell, - copy_data.dofs_per_cell)); - copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell)); - fe_values.reinit (cell, face); - - if (fe_is_system) - // FE has several components - { - boundary_functions.find(cell->face(face)->boundary_id()) - ->second->vector_value_list (fe_values.get_quadrature_points(), - rhs_values_system); - - if (coefficient_is_vector) - // If coefficient is - // vector valued, fill - // all components - coefficient->vector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - else - { - // If a scalar - // function is - // given, update - // the values, if - // not, use the - // default one set - // in the - // constructor above - if (coefficient != 0) - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - // Copy scalar - // values into vector - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - coefficient_vector_values[point] = coefficient_values[point]; - } - - // Special treatment - // for Hdiv and Hcurl - // elements, where only - // the normal or - // tangential component - // should be projected. - std::vector<std::vector<double> > normal_adjustment(fe_values.n_quadrature_points, - std::vector<double>(n_components, 1.)); - - for (unsigned int comp = 0; comp<n_components; ++comp) - { - const FiniteElement<dim,spacedim> &base = fe.base_element(fe.component_to_base_index(comp).first); - const unsigned int bcomp = fe.component_to_base_index(comp).second; - - if (!base.conforms(FiniteElementData<dim>::H1) && - base.conforms(FiniteElementData<dim>::Hdiv)) - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - normal_adjustment[point][comp] = fe_values.normal_vector(point)[bcomp] - * fe_values.normal_vector(point)[bcomp]; - } - - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) - if (fe_is_primitive) - { - for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) - { - if (fe.system_to_component_index(j).first - == fe.system_to_component_index(i).first) - { - copy_data.cell_matrix.back()(i,j) - += weight - * fe_values.shape_value(j,point) - * fe_values.shape_value(i,point) - * coefficient_vector_values[point](fe.system_to_component_index(i).first); - } - } - copy_data.cell_vector.back()(i) += fe_values.shape_value(i,point) - * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first]) - * weight; - } - else - { - for (unsigned int comp=0; comp<n_components; ++comp) - { - for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) - copy_data.cell_matrix.back()(i,j) - += fe_values.shape_value_component(j,point,comp) - * fe_values.shape_value_component(i,point,comp) - * normal_adjustment[point][comp] - * weight * coefficient_vector_values[point](comp); - copy_data.cell_vector.back()(i) += fe_values.shape_value_component(i,point,comp) * - rhs_values_system[point](component_mapping[comp]) - * normal_adjustment[point][comp] - * weight; - } - } - } - } - else - // FE is a scalar one - { - boundary_functions.find(cell->face(face)->boundary_id()) - ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); - - if (coefficient != 0) - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) - { - const double v = fe_values.shape_value(i,point); - for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) - { - const double u = fe_values.shape_value(j,point); - copy_data.cell_matrix.back()(i,j) += (u*v*weight*coefficient_values[point]); - } - copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight; - } - } - } - - - cell->face(face)->get_dof_indices (dofs_on_face_vector); - // for each dof on the cell, have a - // flag whether it is on the face - copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell)); - // check for each of the dofs on this cell - // whether it is on the face - for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) - copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(), - dofs_on_face_vector.end(), - copy_data.dofs[i]) - != - dofs_on_face_vector.end()); - } - } - - template <int dim,int spacedim> - void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim, - spacedim> > const ©_data, - typename FunctionMap<spacedim>::type const &boundary_functions, - std::vector<types::global_dof_index> const &dof_to_boundary_mapping, - SparseMatrix<double> &matrix, - Vector<double> &rhs_vector) - { - // now transfer cell matrix and vector to the whole boundary matrix - // - // in the following: dof[i] holds the global index of the i-th degree of - // freedom on the present cell. If it is also a dof on the boundary, it - // must be a nonzero entry in the dof_to_boundary_mapping and then - // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. - // - // if dof[i] is not on the boundary, it should be zero on the boundary - // therefore on all quadrature points and finally all of its - // entries in the cell matrix and vector should be zero. If not, we - // throw an error (note: because of the evaluation of the shape - // functions only up to machine precision, the term "must be zero" - // really should mean: "should be very small". since this is only an - // assertion and not part of the code, we may choose "very small" - // quite arbitrarily) - // - // the main problem here is that the matrix or vector entry should also - // be zero if the degree of freedom dof[i] is on the boundary, but not - // on the present face, i.e. on another face of the same cell also - // on the boundary. We can therefore not rely on the - // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to - // determine whether dof[i] is a dof on the present face. We do so - // by getting the dofs on the face into @p{dofs_on_face_vector}, - // a vector as always. Usually, searching in a vector is - // inefficient, so we copy the dofs into a set, which enables binary - // searches. - unsigned int pos(0); - for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) - { - // check if this face is on that part of - // the boundary we are interested in - if (boundary_functions.find(copy_data.cell->face(face)->boundary_id()) != - boundary_functions.end()) - { - for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) - { - if (copy_data.dof_is_on_face[pos][i] && - dof_to_boundary_mapping[copy_data.dofs[i]] != numbers::invalid_dof_index) - { - for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j) - if (copy_data.dof_is_on_face[pos][j] && - dof_to_boundary_mapping[copy_data.dofs[j]] != numbers::invalid_dof_index) - { - AssertIsFinite(copy_data.cell_matrix[pos](i,j)); - matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]], - dof_to_boundary_mapping[copy_data.dofs[j]], - copy_data.cell_matrix[pos](i,j)); - } - AssertIsFinite(copy_data.cell_vector[pos](i)); - rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += copy_data.cell_vector[pos](i); - } - } - ++pos; - } - } - } - - - template <> - void - create_boundary_mass_matrix_1<1,3> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<1, - 3> > &/*copy_data*/, - Mapping<1,3> const &, - FiniteElement<1,3> const &, - Quadrature<0> const &, - FunctionMap<3>::type const &/*boundary_functions*/, - Function<3> const *const /*coefficient*/, - std::vector<unsigned int> const &/*component_mapping*/) - { - Assert(false,ExcNotImplemented()); - } - } - - - - - template <int dim, int spacedim> - void - create_boundary_mass_matrix (const Mapping<dim, spacedim> &mapping, - const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim-1> &q, - SparseMatrix<double> &matrix, - const typename FunctionMap<spacedim>::type &boundary_functions, - Vector<double> &rhs_vector, - std::vector<types::global_dof_index> &dof_to_boundary_mapping, - const Function<spacedim> *const coefficient, - std::vector<unsigned int> component_mapping) - { - // what would that be in 1d? the - // identity matrix on the boundary - // dofs? - if (dim == 1) - { - Assert (false, ExcNotImplemented()); - return; - } - - const FiniteElement<dim,spacedim> &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); - - Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), - ExcInternalError()); - Assert (matrix.n() == matrix.m(), ExcInternalError()); - Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); - Assert (boundary_functions.size() != 0, ExcInternalError()); - Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), - ExcInternalError()); - Assert (coefficient ==0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcComponentMismatch()); - - if (component_mapping.size() == 0) - { - AssertDimension (n_components, boundary_functions.begin()->second->n_components); - for (unsigned int i=0; i<n_components; ++i) - component_mapping.push_back(i); - } - else - AssertDimension (n_components, component_mapping.size()); - - MatrixCreator::internal::AssemblerBoundary::Scratch scratch; - MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > copy_data; - - WorkStream::run(dof.begin_active(),dof.end(), - static_cast<std_cxx11::function<void (typename DoFHandler<dim,spacedim>::active_cell_iterator - const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > &)> > - (std_cxx11::bind(&create_boundary_mass_matrix_1<dim,spacedim>,std_cxx11::_1,std_cxx11::_2, - std_cxx11::_3, - std_cxx11::cref(mapping),std_cxx11::cref(fe),std_cxx11::cref(q), - std_cxx11::cref(boundary_functions),coefficient, - std_cxx11::cref(component_mapping))), - static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary - ::CopyData<DoFHandler<dim,spacedim> > const &)> > (std_cxx11::bind( - ©_boundary_mass_matrix_1<dim,spacedim>, - std_cxx11::_1, - std_cxx11::cref(boundary_functions), - std_cxx11::cref(dof_to_boundary_mapping), - std_cxx11::ref(matrix), - std_cxx11::ref(rhs_vector))), - scratch, - copy_data); - } - - - - namespace - { - - template <int dim, int spacedim> - void - create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator const - &cell, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary - ::CopyData<hp::DoFHandler<dim,spacedim> > ©_data, - hp::MappingCollection<dim,spacedim> const &mapping, - hp::FECollection<dim,spacedim> const &fe_collection, - hp::QCollection<dim-1> const &q, - const typename FunctionMap<spacedim>::type &boundary_functions, - Function<spacedim> const *const coefficient, - std::vector<unsigned int> const &component_mapping) - { - const unsigned int n_components = fe_collection.n_components(); - const unsigned int n_function_components = boundary_functions.begin()->second->n_components; - const bool fe_is_system = (n_components != 1); - const FiniteElement<dim,spacedim> &fe = cell->get_fe(); - const unsigned int dofs_per_face = fe.dofs_per_face; - - copy_data.cell = cell; - copy_data.dofs_per_cell = fe.dofs_per_cell; - copy_data.dofs.resize(copy_data.dofs_per_cell); - cell->get_dof_indices (copy_data.dofs); - - - UpdateFlags update_flags = UpdateFlags (update_values | - update_JxW_values | - update_quadrature_points); - hp::FEFaceValues<dim,spacedim> x_fe_values (mapping, fe_collection, q, update_flags); - - // two variables for the coefficient, - // one for the two cases indicated in - // the name - std::vector<double> coefficient_values; - std::vector<Vector<double> > coefficient_vector_values; - - std::vector<double> rhs_values_scalar; - std::vector<Vector<double> > rhs_values_system; - - std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face); - - copy_data.dofs.resize(copy_data.dofs_per_cell); - cell->get_dof_indices (copy_data.dofs); - - // Because CopyData objects are reused and that push_back is used, - // dof_is_on_face, cell_matrix, and cell_vector must be cleared before - // they are reused - copy_data.dof_is_on_face.clear(); - copy_data.cell_matrix.clear(); - copy_data.cell_vector.clear(); - - - for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) - // check if this face is on that part of - // the boundary we are interested in - if (boundary_functions.find(cell->face(face)->boundary_id()) != - boundary_functions.end()) - { - x_fe_values.reinit (cell, face); - - const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values (); - - copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell, - copy_data.dofs_per_cell)); - copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell)); - - if (fe_is_system) - // FE has several components - { - rhs_values_system.resize (fe_values.n_quadrature_points, - Vector<double>(n_function_components)); - boundary_functions.find(cell->face(face)->boundary_id()) - ->second->vector_value_list (fe_values.get_quadrature_points(), - rhs_values_system); - - if (coefficient != 0) - { - if (coefficient->n_components==1) - { - coefficient_values.resize (fe_values.n_quadrature_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) - { - const double v = fe_values.shape_value(i,point); - for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) - if (fe.system_to_component_index(i).first == - fe.system_to_component_index(j).first) - { - const double u = fe_values.shape_value(j,point); - copy_data.cell_matrix.back()(i,j) - += (u * v * weight * coefficient_values[point]); - } - - copy_data.cell_vector.back()(i) += v * - rhs_values_system[point]( - component_mapping[fe.system_to_component_index(i).first]) * weight; - } - } - } - else - { - coefficient_vector_values.resize (fe_values.n_quadrature_points, - Vector<double>(n_components)); - coefficient->vector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) - { - const double v = fe_values.shape_value(i,point); - const unsigned int component_i= - fe.system_to_component_index(i).first; - for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) - if (fe.system_to_component_index(j).first == - component_i) - { - const double u = fe_values.shape_value(j,point); - copy_data.cell_matrix.back()(i,j) += - (u * v * weight * coefficient_vector_values[point](component_i)); - } - copy_data.cell_vector.back()(i) += v * - rhs_values_system[point](component_mapping[component_i]) * weight; - } - } - } - } - else // if (coefficient == 0) - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) - { - const double v = fe_values.shape_value(i,point); - for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) - if (fe.system_to_component_index(i).first == - fe.system_to_component_index(j).first) - { - const double u = fe_values.shape_value(j,point); - copy_data.cell_matrix.back()(i,j) += (u * v * weight); - } - copy_data.cell_vector.back()(i) += v * - rhs_values_system[point]( - fe.system_to_component_index(i).first) * - weight; - } - } - } - else - // FE is a scalar one - { - rhs_values_scalar.resize (fe_values.n_quadrature_points); - boundary_functions.find(cell->face(face)->boundary_id()) - ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); - - if (coefficient != 0) - { - coefficient_values.resize (fe_values.n_quadrature_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) - { - const double v = fe_values.shape_value(i,point); - for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) - { - const double u = fe_values.shape_value(j,point); - copy_data.cell_matrix.back()(i,j) += (u * v * weight * - coefficient_values[point]); - } - copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight; - } - } - } - else - for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) - { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i) - { - const double v = fe_values.shape_value(i,point); - for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j) - { - const double u = fe_values.shape_value(j,point); - copy_data.cell_matrix.back()(i,j) += (u * v * weight); - } - copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] * weight; - } - } - } - - cell->face(face)->get_dof_indices (dofs_on_face_vector, - cell->active_fe_index()); - // for each dof on the cell, have a - // flag whether it is on the face - copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell)); - // check for each of the dofs on this cell - // whether it is on the face - for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) - copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(), - dofs_on_face_vector.end(), - copy_data.dofs[i]) - != - dofs_on_face_vector.end()); - } - } - - - - template <int dim,int spacedim> - void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary - ::CopyData<hp::DoFHandler<dim,spacedim> > const ©_data, - typename FunctionMap<spacedim>::type const &boundary_functions, - std::vector<types::global_dof_index> const &dof_to_boundary_mapping, - SparseMatrix<double> &matrix, - Vector<double> &rhs_vector) - { - // now transfer cell matrix and vector to the whole boundary matrix - // - // in the following: dof[i] holds the global index of the i-th degree of - // freedom on the present cell. If it is also a dof on the boundary, it - // must be a nonzero entry in the dof_to_boundary_mapping and then - // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. - // - // if dof[i] is not on the boundary, it should be zero on the boundary - // therefore on all quadrature points and finally all of its - // entries in the cell matrix and vector should be zero. If not, we - // throw an error (note: because of the evaluation of the shape - // functions only up to machine precision, the term "must be zero" - // really should mean: "should be very small". since this is only an - // assertion and not part of the code, we may choose "very small" - // quite arbitrarily) - // - // the main problem here is that the matrix or vector entry should also - // be zero if the degree of freedom dof[i] is on the boundary, but not - // on the present face, i.e. on another face of the same cell also - // on the boundary. We can therefore not rely on the - // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to - // determine whether dof[i] is a dof on the present face. We do so - // by getting the dofs on the face into @p{dofs_on_face_vector}, - // a vector as always. Usually, searching in a vector is - // inefficient, so we copy the dofs into a set, which enables binary - // searches. - unsigned int pos(0); - for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) - { - // check if this face is on that part of - // the boundary we are interested in - if (boundary_functions.find(copy_data.cell->face(face)->boundary_id()) != - boundary_functions.end()) - { -#ifdef DEBUG - // in debug mode: compute an element in the matrix which is - // guaranteed to belong to a boundary dof. We do this to check that the - // entries in the cell matrix are guaranteed to be zero if the - // respective dof is not on the boundary. Since because of - // round-off, the actual value of the matrix entry may be - // only close to zero, we assert that it is small relative to an element - // which is guaranteed to be nonzero. (absolute smallness does not - // suffice since the size of the domain scales in here) - // - // for this purpose we seek the diagonal of the matrix, where there - // must be an element belonging to the boundary. we take the maximum - // diagonal entry. - types::global_dof_index max_element = static_cast<types::global_dof_index>(0); - for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - Assert (max_element == matrix.n()-1, ExcInternalError()); - - double max_diag_entry = 0; - for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) - if (std::fabs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry) - max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i)); -#endif - - for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i) - for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j) - { - if (copy_data.dof_is_on_face[pos][i] && copy_data.dof_is_on_face[pos][j]) - matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]], - dof_to_boundary_mapping[copy_data.dofs[j]], - copy_data.cell_matrix[pos](i,j)); - else - { - // assume that all shape functions that are nonzero on the boundary - // are also listed in the @p{dof_to_boundary} mapping. if that - // is not the case, then the boundary mass matrix does not - // make that much sense anyway, as it only contains entries for - // parts of the functions living on the boundary - // - // these, we may compare here for relative smallness of all - // entries in the local matrix which are not taken over to - // the global one - Assert (std::fabs(copy_data.cell_matrix[pos](i,j)) <= 1e-10 * max_diag_entry, - ExcInternalError ()); - } - } - - for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j) - if (copy_data.dof_is_on_face[pos][j]) - rhs_vector(dof_to_boundary_mapping[copy_data.dofs[j]]) += copy_data.cell_vector[pos](j); - else - { - // compare here for relative - // smallness - Assert (std::fabs(copy_data.cell_vector[pos](j)) <= 1e-10 * max_diag_entry, - ExcInternalError()); - } - ++pos; - } - } - } - } - - - - template <int dim, int spacedim> - void create_boundary_mass_matrix (const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim-1> &q, - SparseMatrix<double> &matrix, - const typename FunctionMap<spacedim>::type &rhs, - Vector<double> &rhs_vector, - std::vector<types::global_dof_index> &dof_to_boundary_mapping, - const Function<spacedim> *const a, - std::vector<unsigned int> component_mapping) - { - create_boundary_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q, - matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); - } - - - - template <int dim, int spacedim> - void - create_boundary_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping, - const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim-1> &q, - SparseMatrix<double> &matrix, - const typename FunctionMap<spacedim>::type &boundary_functions, - Vector<double> &rhs_vector, - std::vector<types::global_dof_index> &dof_to_boundary_mapping, - const Function<spacedim> *const coefficient, - std::vector<unsigned int> component_mapping) - { - // what would that be in 1d? the - // identity matrix on the boundary - // dofs? - if (dim == 1) - { - Assert (false, ExcNotImplemented()); - return; - } - - const hp::FECollection<dim,spacedim> &fe_collection = dof.get_fe(); - const unsigned int n_components = fe_collection.n_components(); - - Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), - ExcInternalError()); - Assert (matrix.n() == matrix.m(), ExcInternalError()); - Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); - Assert (boundary_functions.size() != 0, ExcInternalError()); - Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), - ExcInternalError()); - Assert (coefficient ==0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcComponentMismatch()); - - if (component_mapping.size() == 0) - { - AssertDimension (n_components, boundary_functions.begin()->second->n_components); - for (unsigned int i=0; i<n_components; ++i) - component_mapping.push_back(i); - } - else - AssertDimension (n_components, component_mapping.size()); - - MatrixCreator::internal::AssemblerBoundary::Scratch scratch; - MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > copy_data; - - WorkStream::run(dof.begin_active(),dof.end(), - static_cast<std_cxx11::function<void (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator - const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > &)> > - (std_cxx11::bind( &create_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx11::_1,std_cxx11::_2, - std_cxx11::_3, - std_cxx11::cref(mapping),std_cxx11::cref(fe_collection),std_cxx11::cref(q), - std_cxx11::cref(boundary_functions),coefficient, - std_cxx11::cref(component_mapping))), - static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary - ::CopyData<hp::DoFHandler<dim,spacedim> > const &)> > ( - std_cxx11::bind( ©_hp_boundary_mass_matrix_1<dim,spacedim>, - std_cxx11::_1, - std_cxx11::cref(boundary_functions), - std_cxx11::cref(dof_to_boundary_mapping), - std_cxx11::ref(matrix), - std_cxx11::ref(rhs_vector))), - scratch, - copy_data); - } - - - - - template <int dim, int spacedim> - void create_boundary_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim-1> &q, - SparseMatrix<double> &matrix, - const typename FunctionMap<spacedim>::type &rhs, - Vector<double> &rhs_vector, - std::vector<types::global_dof_index> &dof_to_boundary_mapping, - const Function<spacedim> *const a, - std::vector<unsigned int> component_mapping) - { - create_boundary_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, - matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); - } - - - - template <int dim, int spacedim> - void create_laplace_matrix (const Mapping<dim, spacedim> &mapping, - const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim> &q, - SparseMatrix<double> &matrix, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - hp::FECollection<dim,spacedim> fe_collection (dof.get_fe()); - hp::QCollection<dim> q_collection (q); - hp::MappingCollection<dim,spacedim> mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> - assembler_data (fe_collection, - update_gradients | update_JxW_values | - (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), - coefficient, /*rhs_function=*/0, - q_collection, mapping_collection); - MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; - - WorkStream::run (dof.begin_active(), - static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global<double,SparseMatrix<double>, Vector<double> >, - std_cxx11::_1, - &matrix, - (Vector<double> *)NULL), - assembler_data, - copy_data); - } - - - - template <int dim, int spacedim> - void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim> &q, - SparseMatrix<double> &matrix, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, - dof, q, matrix, coefficient, constraints); - } - - - - template <int dim, int spacedim> - void create_laplace_matrix (const Mapping<dim, spacedim> &mapping, - const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim> &q, - SparseMatrix<double> &matrix, - const Function<spacedim> &rhs, - Vector<double> &rhs_vector, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - hp::FECollection<dim,spacedim> fe_collection (dof.get_fe()); - hp::QCollection<dim> q_collection (q); - hp::MappingCollection<dim,spacedim> mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> - assembler_data (fe_collection, - update_gradients | update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q_collection, mapping_collection); - MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; - - WorkStream::run (dof.begin_active(), - static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global<double,SparseMatrix<double>, Vector<double> >, - std_cxx11::_1, - &matrix, - &rhs_vector), - assembler_data, - copy_data); - } - - - - template <int dim, int spacedim> - void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof, - const Quadrature<dim> &q, - SparseMatrix<double> &matrix, - const Function<spacedim> &rhs, - Vector<double> &rhs_vector, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q, - matrix, rhs, rhs_vector, coefficient, constraints); - } - - - - template <int dim, int spacedim> - void create_laplace_matrix (const hp::MappingCollection<dim,spacedim> &mapping, - const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim> &q, - SparseMatrix<double> &matrix, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> - assembler_data (dof.get_fe(), - update_gradients | update_JxW_values | - (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), - coefficient, /*rhs_function=*/0, - q, mapping); - MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; - - WorkStream::run (dof.begin_active(), - static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global<double,SparseMatrix<double>, Vector<double> >, - std_cxx11::_1, - &matrix, - (Vector<double> *)0), - assembler_data, - copy_data); - } - - - - template <int dim, int spacedim> - void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim> &q, - SparseMatrix<double> &matrix, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, - matrix, coefficient, constraints); - } - - - - template <int dim, int spacedim> - void create_laplace_matrix (const hp::MappingCollection<dim,spacedim> &mapping, - const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim> &q, - SparseMatrix<double> &matrix, - const Function<spacedim> &rhs, - Vector<double> &rhs_vector, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim> - assembler_data (dof.get_fe(), - update_gradients | update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q, mapping); - MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; - - WorkStream::run (dof.begin_active(), - static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global<double,SparseMatrix<double>, Vector<double> >, - std_cxx11::_1, - &matrix, - &rhs_vector), - assembler_data, - copy_data); - } - - - - template <int dim, int spacedim> - void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof, - const hp::QCollection<dim> &q, - SparseMatrix<double> &matrix, - const Function<spacedim> &rhs, - Vector<double> &rhs_vector, - const Function<spacedim> *const coefficient, - const ConstraintMatrix &constraints) - { - create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, - matrix, rhs, rhs_vector, coefficient, constraints); - } - -} // namespace MatrixCreator - - namespace MatrixTools { namespace @@ -2316,462 +489,6 @@ namespace MatrixTools -#ifdef DEAL_II_WITH_PETSC - - namespace internal - { - namespace PETScWrappers - { - template <typename PETScMatrix, typename PETScVector> - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - PETScMatrix &matrix, - PETScVector &solution, - PETScVector &right_hand_side, - const bool eliminate_columns) - { - (void)eliminate_columns; - Assert (eliminate_columns == false, ExcNotImplemented()); - - Assert (matrix.n() == right_hand_side.size(), - ExcDimensionMismatch(matrix.n(), right_hand_side.size())); - Assert (matrix.n() == solution.size(), - ExcDimensionMismatch(matrix.n(), solution.size())); - - // if no boundary values are to be applied, then - // jump straight to the compress() calls that we still have - // to perform because they are collective operations - if (boundary_values.size() > 0) - { - const std::pair<types::global_dof_index, types::global_dof_index> local_range - = matrix.local_range(); - Assert (local_range == right_hand_side.local_range(), - ExcInternalError()); - Assert (local_range == solution.local_range(), - ExcInternalError()); - - // determine the first nonzero diagonal - // entry from within the part of the - // matrix that we can see. if we can't - // find such an entry, take one - PetscScalar average_nonzero_diagonal_entry = 1; - for (types::global_dof_index i=local_range.first; i<local_range.second; ++i) - if (matrix.diag_element(i) != PetscScalar ()) - { - average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i)); - break; - } - - // figure out which rows of the matrix we - // have to eliminate on this processor - std::vector<types::global_dof_index> constrained_rows; - for (std::map<types::global_dof_index,double>::const_iterator - dof = boundary_values.begin(); - dof != boundary_values.end(); - ++dof) - if ((dof->first >= local_range.first) && - (dof->first < local_range.second)) - constrained_rows.push_back (dof->first); - - // then eliminate these rows and set - // their diagonal entry to what we have - // determined above. note that for petsc - // matrices interleaving read with write - // operations is very expensive. thus, we - // here always replace the diagonal - // element, rather than first checking - // whether it is nonzero and in that case - // preserving it. this is different from - // the case of deal.II sparse matrices - // treated in the other functions. - matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry); - - std::vector<types::global_dof_index> indices; - std::vector<PetscScalar> solution_values; - for (std::map<types::global_dof_index,double>::const_iterator - dof = boundary_values.begin(); - dof != boundary_values.end(); - ++dof) - if ((dof->first >= local_range.first) && - (dof->first < local_range.second)) - { - indices.push_back (dof->first); - solution_values.push_back (dof->second); - } - solution.set (indices, solution_values); - - // now also set appropriate values for - // the rhs - for (unsigned int i=0; i<solution_values.size(); ++i) - solution_values[i] *= average_nonzero_diagonal_entry; - - right_hand_side.set (indices, solution_values); - } - else - { - // clear_rows() is a collective operation so we still have to call - // it: - std::vector<types::global_dof_index> constrained_rows; - matrix.clear_rows (constrained_rows, 1.); - } - - // clean up - solution.compress (VectorOperation::insert); - right_hand_side.compress (VectorOperation::insert); - } - } - } - - - - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - PETScWrappers::SparseMatrix &matrix, - PETScWrappers::Vector &solution, - PETScWrappers::Vector &right_hand_side, - const bool eliminate_columns) - { - // simply redirect to the generic function - // used for both petsc matrix types - internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution, - right_hand_side, eliminate_columns); - } - - - - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - PETScWrappers::MPI::SparseMatrix &matrix, - PETScWrappers::MPI::Vector &solution, - PETScWrappers::MPI::Vector &right_hand_side, - const bool eliminate_columns) - { - // simply redirect to the generic function - // used for both petsc matrix types - internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution, - right_hand_side, eliminate_columns); - } - - - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - PETScWrappers::MPI::BlockSparseMatrix &matrix, - PETScWrappers::MPI::BlockVector &solution, - PETScWrappers::MPI::BlockVector &right_hand_side, - const bool eliminate_columns) - { - Assert (matrix.n() == right_hand_side.size(), - ExcDimensionMismatch(matrix.n(), right_hand_side.size())); - Assert (matrix.n() == solution.size(), - ExcDimensionMismatch(matrix.n(), solution.size())); - Assert (matrix.n_block_rows() == matrix.n_block_cols(), - ExcNotQuadratic()); - - const unsigned int n_blocks = matrix.n_block_rows(); - - // We need to find the subdivision - // into blocks for the boundary values. - // To this end, generate a vector of - // maps with the respective indices. - std::vector<std::map<dealii::types::global_dof_index,double> > block_boundary_values(n_blocks); - { - int block = 0; - dealii::types::global_dof_index offset = 0; - for (std::map<types::global_dof_index,double>::const_iterator - dof = boundary_values.begin(); - dof != boundary_values.end(); - ++dof) - { - if (dof->first >= matrix.block(block,0).m() + offset) - { - offset += matrix.block(block,0).m(); - block++; - } - const types::global_dof_index index = dof->first - offset; - block_boundary_values[block].insert(std::pair<types::global_dof_index, double> (index,dof->second)); - } - } - - // Now call the non-block variants on - // the diagonal subblocks and the - // solution/rhs. - for (unsigned int block=0; block<n_blocks; ++block) - internal::PETScWrappers::apply_boundary_values(block_boundary_values[block], - matrix.block(block,block), - solution.block(block), - right_hand_side.block(block), - eliminate_columns); - - // Finally, we need to do something - // about the off-diagonal matrices. This - // is luckily not difficult. Just clear - // the whole row. - for (unsigned int block_m=0; block_m<n_blocks; ++block_m) - { - const std::pair<types::global_dof_index, types::global_dof_index> local_range - = matrix.block(block_m,0).local_range(); - - std::vector<types::global_dof_index> constrained_rows; - for (std::map<types::global_dof_index,double>::const_iterator - dof = block_boundary_values[block_m].begin(); - dof != block_boundary_values[block_m].end(); - ++dof) - if ((dof->first >= local_range.first) && - (dof->first < local_range.second)) - constrained_rows.push_back (dof->first); - - for (unsigned int block_n=0; block_n<n_blocks; ++block_n) - if (block_m != block_n) - matrix.block(block_m,block_n).clear_rows(constrained_rows); - } - } - -#endif - - - -#ifdef DEAL_II_WITH_TRILINOS - - namespace internal - { - namespace TrilinosWrappers - { - template <typename TrilinosMatrix, typename TrilinosVector> - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - TrilinosMatrix &matrix, - TrilinosVector &solution, - TrilinosVector &right_hand_side, - const bool eliminate_columns) - { - Assert (eliminate_columns == false, ExcNotImplemented()); - (void)eliminate_columns; - - Assert (matrix.n() == right_hand_side.size(), - ExcDimensionMismatch(matrix.n(), right_hand_side.size())); - Assert (matrix.n() == solution.size(), - ExcDimensionMismatch(matrix.m(), solution.size())); - - // if no boundary values are to be applied, then - // jump straight to the compress() calls that we still have - // to perform because they are collective operations - if (boundary_values.size() > 0) - { - const std::pair<types::global_dof_index, types::global_dof_index> local_range - = matrix.local_range(); - Assert (local_range == right_hand_side.local_range(), - ExcInternalError()); - Assert (local_range == solution.local_range(), - ExcInternalError()); - - // determine the first nonzero diagonal - // entry from within the part of the - // matrix that we can see. if we can't - // find such an entry, take one - TrilinosScalar average_nonzero_diagonal_entry = 1; - for (types::global_dof_index i=local_range.first; i<local_range.second; ++i) - if (matrix.diag_element(i) != 0) - { - average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i)); - break; - } - - // figure out which rows of the matrix we - // have to eliminate on this processor - std::vector<types::global_dof_index> constrained_rows; - for (std::map<types::global_dof_index,double>::const_iterator - dof = boundary_values.begin(); - dof != boundary_values.end(); - ++dof) - if ((dof->first >= local_range.first) && - (dof->first < local_range.second)) - constrained_rows.push_back (dof->first); - - // then eliminate these rows and - // set their diagonal entry to - // what we have determined - // above. if the value already is - // nonzero, it will be preserved, - // in accordance with the basic - // matrix classes in deal.II. - matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry); - - std::vector<types::global_dof_index> indices; - std::vector<TrilinosScalar> solution_values; - for (std::map<types::global_dof_index,double>::const_iterator - dof = boundary_values.begin(); - dof != boundary_values.end(); - ++dof) - if ((dof->first >= local_range.first) && - (dof->first < local_range.second)) - { - indices.push_back (dof->first); - solution_values.push_back (dof->second); - } - solution.set (indices, solution_values); - - // now also set appropriate - // values for the rhs - for (unsigned int i=0; i<solution_values.size(); ++i) - solution_values[i] *= matrix.diag_element(indices[i]); - - right_hand_side.set (indices, solution_values); - } - else - { - // clear_rows() is a collective operation so we still have to call - // it: - std::vector<types::global_dof_index> constrained_rows; - matrix.clear_rows (constrained_rows, 1.); - } - - // clean up - matrix.compress (VectorOperation::insert); - solution.compress (VectorOperation::insert); - right_hand_side.compress (VectorOperation::insert); - } - - - - template <typename TrilinosMatrix, typename TrilinosBlockVector> - void - apply_block_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - TrilinosMatrix &matrix, - TrilinosBlockVector &solution, - TrilinosBlockVector &right_hand_side, - const bool eliminate_columns) - { - Assert (eliminate_columns == false, ExcNotImplemented()); - - Assert (matrix.n() == right_hand_side.size(), - ExcDimensionMismatch(matrix.n(), right_hand_side.size())); - Assert (matrix.n() == solution.size(), - ExcDimensionMismatch(matrix.n(), solution.size())); - Assert (matrix.n_block_rows() == matrix.n_block_cols(), - ExcNotQuadratic()); - - const unsigned int n_blocks = matrix.n_block_rows(); - - // We need to find the subdivision - // into blocks for the boundary values. - // To this end, generate a vector of - // maps with the respective indices. - std::vector<std::map<types::global_dof_index,double> > block_boundary_values(n_blocks); - { - int block=0; - types::global_dof_index offset = 0; - for (std::map<types::global_dof_index,double>::const_iterator - dof = boundary_values.begin(); - dof != boundary_values.end(); - ++dof) - { - if (dof->first >= matrix.block(block,0).m() + offset) - { - offset += matrix.block(block,0).m(); - block++; - } - const types::global_dof_index index = dof->first - offset; - block_boundary_values[block].insert( - std::pair<types::global_dof_index, double> (index,dof->second)); - } - } - - // Now call the non-block variants on - // the diagonal subblocks and the - // solution/rhs. - for (unsigned int block=0; block<n_blocks; ++block) - TrilinosWrappers::apply_boundary_values(block_boundary_values[block], - matrix.block(block,block), - solution.block(block), - right_hand_side.block(block), - eliminate_columns); - - // Finally, we need to do something - // about the off-diagonal matrices. This - // is luckily not difficult. Just clear - // the whole row. - for (unsigned int block_m=0; block_m<n_blocks; ++block_m) - { - const std::pair<types::global_dof_index, types::global_dof_index> local_range - = matrix.block(block_m,0).local_range(); - - std::vector<types::global_dof_index> constrained_rows; - for (std::map<types::global_dof_index,double>::const_iterator - dof = block_boundary_values[block_m].begin(); - dof != block_boundary_values[block_m].end(); - ++dof) - if ((dof->first >= local_range.first) && - (dof->first < local_range.second)) - constrained_rows.push_back (dof->first); - - for (unsigned int block_n=0; block_n<n_blocks; ++block_n) - if (block_m != block_n) - matrix.block(block_m,block_n).clear_rows(constrained_rows); - } - } - } - } - - - - - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - TrilinosWrappers::SparseMatrix &matrix, - TrilinosWrappers::Vector &solution, - TrilinosWrappers::Vector &right_hand_side, - const bool eliminate_columns) - { - // simply redirect to the generic function - // used for both trilinos matrix types - internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution, - right_hand_side, eliminate_columns); - } - - - - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - TrilinosWrappers::SparseMatrix &matrix, - TrilinosWrappers::MPI::Vector &solution, - TrilinosWrappers::MPI::Vector &right_hand_side, - const bool eliminate_columns) - { - // simply redirect to the generic function - // used for both trilinos matrix types - internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution, - right_hand_side, eliminate_columns); - } - - - - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - TrilinosWrappers::BlockSparseMatrix &matrix, - TrilinosWrappers::BlockVector &solution, - TrilinosWrappers::BlockVector &right_hand_side, - const bool eliminate_columns) - { - internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix, - solution, right_hand_side, - eliminate_columns); - } - - - - void - apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, - TrilinosWrappers::BlockSparseMatrix &matrix, - TrilinosWrappers::MPI::BlockVector &solution, - TrilinosWrappers::MPI::BlockVector &right_hand_side, - const bool eliminate_columns) - { - internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix, - solution, right_hand_side, - eliminate_columns); - } - -#endif @@ -2894,37 +611,5 @@ namespace MatrixTools // explicit instantiations #include "matrix_tools.inst" -namespace MatrixTools -{ - template - void - apply_boundary_values<double> (const std::map<types::global_dof_index,double> &boundary_values, - SparseMatrix<double> &matrix, - Vector<double> &solution, - Vector<double> &right_hand_side, - const bool eliminate_columns); - template - void - apply_boundary_values<float> (const std::map<types::global_dof_index,double> &boundary_values, - SparseMatrix<float> &matrix, - Vector<float> &solution, - Vector<float> &right_hand_side, - const bool eliminate_columns); - - template - void - apply_boundary_values<double> (const std::map<types::global_dof_index,double> &boundary_values, - BlockSparseMatrix<double> &matrix, - BlockVector<double> &solution, - BlockVector<double> &right_hand_side, - const bool eliminate_columns); - template - void - apply_boundary_values<float> (const std::map<types::global_dof_index,double> &boundary_values, - BlockSparseMatrix<float> &matrix, - BlockVector<float> &solution, - BlockVector<float> &right_hand_side, - const bool eliminate_columns); -} DEAL_II_NAMESPACE_CLOSE diff --git a/source/numerics/matrix_tools.inst.in b/source/numerics/matrix_tools.inst.in index a2e9346dc5..3516f400f3 100644 --- a/source/numerics/matrix_tools.inst.in +++ b/source/numerics/matrix_tools.inst.in @@ -13,306 +13,23 @@ // // --------------------------------------------------------------------- -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) - { -#if deal_II_dimension <= deal_II_space_dimension - -// non-hp version of create_mass_matrix - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> - (const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping, - const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> - (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> - (const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping, - const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_space_dimension> &rhs, - Vector<double> &rhs_vector, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> - (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_space_dimension> &rhs, - Vector<double> &rhs_vector, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - - - template - void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension> - (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension-1> &q, - SparseMatrix<double> &matrix, - const FunctionMap<deal_II_space_dimension>::type &rhs, - Vector<double> &rhs_vector, - std::vector<types::global_dof_index> &dof_to_boundary_mapping, - const Function<deal_II_space_dimension> * const a, - std::vector<unsigned int>); - - template - void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension> - (const Mapping<deal_II_dimension,deal_II_space_dimension> &, - const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension-1> &q, - SparseMatrix<double> &matrix, - const FunctionMap<deal_II_space_dimension>::type &rhs, - Vector<double> &rhs_vector, - std::vector<types::global_dof_index> &dof_to_boundary_mapping, - const Function<deal_II_space_dimension> * const a, - std::vector<unsigned int>); - - template - void - MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension> - (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>&, - const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&, - const hp::QCollection<deal_II_dimension-1>&, - SparseMatrix<double>&, - const FunctionMap<deal_II_space_dimension>::type&, - Vector<double>&, - std::vector<types::global_dof_index>&, - const Function<deal_II_space_dimension> * const, - std::vector<unsigned int>); - - template - void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension> - (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&, - const hp::QCollection<deal_II_dimension-1>&, - SparseMatrix<double>&, - const FunctionMap<deal_II_space_dimension>::type&, - Vector<double>&, - std::vector<types::global_dof_index>&, - const Function<deal_II_space_dimension> * const, - std::vector<unsigned int>); - -// same for float - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> - (const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping, - const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<float> &matrix, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> - (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<float> &matrix, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> - (const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping, - const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<float> &matrix, - const Function<deal_II_space_dimension> &rhs, - Vector<float> &rhs_vector, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> - (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<float> &matrix, - const Function<deal_II_space_dimension> &rhs, - Vector<float> &rhs_vector, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - -#endif - } - -//TODO[SP]: replace <deal_II_dimension> by <deal_II_dimension, deal_II_space_dimension> -// where applicable and move to codimension cases above also when applicable -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +for (number: REAL_SCALARS) { -// hp versions of functions -#if deal_II_dimension <= deal_II_space_dimension - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> - (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping, - const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension> - (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping, - const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_space_dimension> &rhs, - Vector<double> &rhs_vector, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - - -#endif - -#if deal_II_dimension == deal_II_space_dimension + void MatrixTools::apply_boundary_values + (const std::map<types::global_dof_index,double> &boundary_values, + SparseMatrix<number> &matrix, + Vector<number> &solution, + Vector<number> &right_hand_side, + const bool eliminate_columns); - template - void MatrixCreator::create_mass_matrix<deal_II_dimension> - (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - - template - void MatrixCreator::create_mass_matrix<deal_II_dimension> - (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_space_dimension> &rhs, - Vector<double> &rhs_vector, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - -// non-hp versions of create_laplace_matrix - template - void MatrixCreator::create_laplace_matrix<deal_II_dimension> - (const DoFHandler<deal_II_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix<deal_II_dimension> - (const Mapping<deal_II_dimension> &mapping, - const DoFHandler<deal_II_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix<deal_II_dimension> - (const Mapping<deal_II_dimension> &mapping, - const DoFHandler<deal_II_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_dimension> &rhs, - Vector<double> &rhs_vector, - const Function<deal_II_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix<deal_II_dimension> - (const DoFHandler<deal_II_dimension> &dof, - const Quadrature<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_dimension> &rhs, - Vector<double> &rhs_vector, - const Function<deal_II_dimension> * const coefficient, - const ConstraintMatrix &constraints); - -// hp versions of create_laplace_matrix - template - void MatrixCreator::create_laplace_matrix<deal_II_dimension> - (const hp::DoFHandler<deal_II_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix<deal_II_dimension> - (const hp::MappingCollection<deal_II_dimension> &mapping, - const hp::DoFHandler<deal_II_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix<deal_II_dimension> - (const hp::MappingCollection<deal_II_dimension> &mapping, - const hp::DoFHandler<deal_II_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_dimension> &rhs, - Vector<double> &rhs_vector, - const Function<deal_II_dimension> * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix<deal_II_dimension> - (const hp::DoFHandler<deal_II_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<double> &matrix, - const Function<deal_II_dimension> &rhs, - Vector<double> &rhs_vector, - const Function<deal_II_dimension> * const coefficient, - const ConstraintMatrix &constraints); - -#endif - -// same for float -#if deal_II_dimension <= deal_II_space_dimension - - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> - (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping, - const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<float> &matrix, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - - template - void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension> - (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension> &mapping, - const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<float> &matrix, - const Function<deal_II_space_dimension> &rhs, - Vector<float> &rhs_vector, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - - -#endif - -#if deal_II_dimension == deal_II_space_dimension - - - template - void MatrixCreator::create_mass_matrix<deal_II_dimension> - (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<float> &matrix, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - - template - void MatrixCreator::create_mass_matrix<deal_II_dimension> - (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof, - const hp::QCollection<deal_II_dimension> &q, - SparseMatrix<float> &matrix, - const Function<deal_II_space_dimension> &rhs, - Vector<float> &rhs_vector, - const Function<deal_II_space_dimension> * const coefficient, - const ConstraintMatrix &constraints); - -#endif + void MatrixTools::apply_boundary_values + (const std::map<types::global_dof_index,double> &boundary_values, + BlockSparseMatrix<number> &matrix, + BlockVector<number> &solution, + BlockVector<number> &right_hand_side, + const bool eliminate_columns); } - diff --git a/source/numerics/matrix_tools_once.cc b/source/numerics/matrix_tools_once.cc new file mode 100644 index 0000000000..d9bd2d1f18 --- /dev/null +++ b/source/numerics/matrix_tools_once.cc @@ -0,0 +1,520 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2015 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. +// +// --------------------------------------------------------------------- + +#include <deal.II/base/function.h> +#include <deal.II/base/quadrature.h> +#include <deal.II/base/work_stream.h> +#include <deal.II/base/geometry_info.h> +#include <deal.II/base/quadrature.h> +#include <deal.II/dofs/dof_handler.h> +#include <deal.II/dofs/dof_accessor.h> +#include <deal.II/dofs/dof_tools.h> +#include <deal.II/fe/fe.h> +#include <deal.II/fe/fe_values.h> +#include <deal.II/fe/mapping_q1.h> +#include <deal.II/grid/tria_iterator.h> +#include <deal.II/hp/fe_values.h> +#include <deal.II/hp/mapping_collection.h> +#include <deal.II/numerics/matrix_tools.h> +#include <deal.II/lac/vector.h> +#include <deal.II/lac/block_vector.h> +#include <deal.II/lac/full_matrix.h> +#include <deal.II/lac/sparse_matrix.h> +#include <deal.II/lac/block_sparse_matrix.h> + +#ifdef DEAL_II_WITH_PETSC +# include <deal.II/lac/petsc_parallel_sparse_matrix.h> +# include <deal.II/lac/petsc_sparse_matrix.h> +# include <deal.II/lac/petsc_parallel_vector.h> +# include <deal.II/lac/petsc_vector.h> +# include <deal.II/lac/petsc_parallel_block_sparse_matrix.h> +#endif + +#ifdef DEAL_II_WITH_TRILINOS +# include <deal.II/lac/trilinos_sparse_matrix.h> +# include <deal.II/lac/trilinos_vector.h> +# include <deal.II/lac/trilinos_block_sparse_matrix.h> +# include <deal.II/lac/trilinos_block_vector.h> +#endif + +#include <algorithm> + + +#include <algorithm> +#include <set> +#include <cmath> + + +DEAL_II_NAMESPACE_OPEN + +#ifdef DEAL_II_WITH_PETSC + +namespace internal +{ + namespace PETScWrappers + { + template <typename PETScMatrix, typename PETScVector> + void + apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + PETScMatrix &matrix, + PETScVector &solution, + PETScVector &right_hand_side, + const bool eliminate_columns) + { + (void)eliminate_columns; + Assert (eliminate_columns == false, ExcNotImplemented()); + + Assert (matrix.n() == right_hand_side.size(), + ExcDimensionMismatch(matrix.n(), right_hand_side.size())); + Assert (matrix.n() == solution.size(), + ExcDimensionMismatch(matrix.n(), solution.size())); + + // if no boundary values are to be applied, then + // jump straight to the compress() calls that we still have + // to perform because they are collective operations + if (boundary_values.size() > 0) + { + const std::pair<types::global_dof_index, types::global_dof_index> local_range + = matrix.local_range(); + Assert (local_range == right_hand_side.local_range(), + ExcInternalError()); + Assert (local_range == solution.local_range(), + ExcInternalError()); + + // determine the first nonzero diagonal + // entry from within the part of the + // matrix that we can see. if we can't + // find such an entry, take one + PetscScalar average_nonzero_diagonal_entry = 1; + for (types::global_dof_index i=local_range.first; i<local_range.second; ++i) + if (matrix.diag_element(i) != PetscScalar ()) + { + average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i)); + break; + } + + // figure out which rows of the matrix we + // have to eliminate on this processor + std::vector<types::global_dof_index> constrained_rows; + for (std::map<types::global_dof_index,double>::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + constrained_rows.push_back (dof->first); + + // then eliminate these rows and set + // their diagonal entry to what we have + // determined above. note that for petsc + // matrices interleaving read with write + // operations is very expensive. thus, we + // here always replace the diagonal + // element, rather than first checking + // whether it is nonzero and in that case + // preserving it. this is different from + // the case of deal.II sparse matrices + // treated in the other functions. + matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry); + + std::vector<types::global_dof_index> indices; + std::vector<PetscScalar> solution_values; + for (std::map<types::global_dof_index,double>::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + { + indices.push_back (dof->first); + solution_values.push_back (dof->second); + } + solution.set (indices, solution_values); + + // now also set appropriate values for + // the rhs + for (unsigned int i=0; i<solution_values.size(); ++i) + solution_values[i] *= average_nonzero_diagonal_entry; + + right_hand_side.set (indices, solution_values); + } + else + { + // clear_rows() is a collective operation so we still have to call + // it: + std::vector<types::global_dof_index> constrained_rows; + matrix.clear_rows (constrained_rows, 1.); + } + + // clean up + solution.compress (VectorOperation::insert); + right_hand_side.compress (VectorOperation::insert); + } + } +} + + + +void +apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + PETScWrappers::SparseMatrix &matrix, + PETScWrappers::Vector &solution, + PETScWrappers::Vector &right_hand_side, + const bool eliminate_columns) +{ + // simply redirect to the generic function + // used for both petsc matrix types + internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution, + right_hand_side, eliminate_columns); +} + + + +void +apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + PETScWrappers::MPI::SparseMatrix &matrix, + PETScWrappers::MPI::Vector &solution, + PETScWrappers::MPI::Vector &right_hand_side, + const bool eliminate_columns) +{ + // simply redirect to the generic function + // used for both petsc matrix types + internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution, + right_hand_side, eliminate_columns); +} + + +void +apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + PETScWrappers::MPI::BlockSparseMatrix &matrix, + PETScWrappers::MPI::BlockVector &solution, + PETScWrappers::MPI::BlockVector &right_hand_side, + const bool eliminate_columns) +{ + Assert (matrix.n() == right_hand_side.size(), + ExcDimensionMismatch(matrix.n(), right_hand_side.size())); + Assert (matrix.n() == solution.size(), + ExcDimensionMismatch(matrix.n(), solution.size())); + Assert (matrix.n_block_rows() == matrix.n_block_cols(), + ExcNotQuadratic()); + + const unsigned int n_blocks = matrix.n_block_rows(); + + // We need to find the subdivision + // into blocks for the boundary values. + // To this end, generate a vector of + // maps with the respective indices. + std::vector<std::map<dealii::types::global_dof_index,double> > block_boundary_values(n_blocks); + { + int block = 0; + dealii::types::global_dof_index offset = 0; + for (std::map<types::global_dof_index,double>::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + { + if (dof->first >= matrix.block(block,0).m() + offset) + { + offset += matrix.block(block,0).m(); + block++; + } + const types::global_dof_index index = dof->first - offset; + block_boundary_values[block].insert(std::pair<types::global_dof_index, double> (index,dof->second)); + } + } + + // Now call the non-block variants on + // the diagonal subblocks and the + // solution/rhs. + for (unsigned int block=0; block<n_blocks; ++block) + internal::PETScWrappers::apply_boundary_values(block_boundary_values[block], + matrix.block(block,block), + solution.block(block), + right_hand_side.block(block), + eliminate_columns); + + // Finally, we need to do something + // about the off-diagonal matrices. This + // is luckily not difficult. Just clear + // the whole row. + for (unsigned int block_m=0; block_m<n_blocks; ++block_m) + { + const std::pair<types::global_dof_index, types::global_dof_index> local_range + = matrix.block(block_m,0).local_range(); + + std::vector<types::global_dof_index> constrained_rows; + for (std::map<types::global_dof_index,double>::const_iterator + dof = block_boundary_values[block_m].begin(); + dof != block_boundary_values[block_m].end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + constrained_rows.push_back (dof->first); + + for (unsigned int block_n=0; block_n<n_blocks; ++block_n) + if (block_m != block_n) + matrix.block(block_m,block_n).clear_rows(constrained_rows); + } +} + +#endif + + + +#ifdef DEAL_II_WITH_TRILINOS + +namespace internal +{ + namespace TrilinosWrappers + { + template <typename TrilinosMatrix, typename TrilinosVector> + void + apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + TrilinosMatrix &matrix, + TrilinosVector &solution, + TrilinosVector &right_hand_side, + const bool eliminate_columns) + { + Assert (eliminate_columns == false, ExcNotImplemented()); + (void)eliminate_columns; + + Assert (matrix.n() == right_hand_side.size(), + ExcDimensionMismatch(matrix.n(), right_hand_side.size())); + Assert (matrix.n() == solution.size(), + ExcDimensionMismatch(matrix.m(), solution.size())); + + // if no boundary values are to be applied, then + // jump straight to the compress() calls that we still have + // to perform because they are collective operations + if (boundary_values.size() > 0) + { + const std::pair<types::global_dof_index, types::global_dof_index> local_range + = matrix.local_range(); + Assert (local_range == right_hand_side.local_range(), + ExcInternalError()); + Assert (local_range == solution.local_range(), + ExcInternalError()); + + // determine the first nonzero diagonal + // entry from within the part of the + // matrix that we can see. if we can't + // find such an entry, take one + TrilinosScalar average_nonzero_diagonal_entry = 1; + for (types::global_dof_index i=local_range.first; i<local_range.second; ++i) + if (matrix.diag_element(i) != 0) + { + average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i)); + break; + } + + // figure out which rows of the matrix we + // have to eliminate on this processor + std::vector<types::global_dof_index> constrained_rows; + for (std::map<types::global_dof_index,double>::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + constrained_rows.push_back (dof->first); + + // then eliminate these rows and + // set their diagonal entry to + // what we have determined + // above. if the value already is + // nonzero, it will be preserved, + // in accordance with the basic + // matrix classes in deal.II. + matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry); + + std::vector<types::global_dof_index> indices; + std::vector<TrilinosScalar> solution_values; + for (std::map<types::global_dof_index,double>::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + { + indices.push_back (dof->first); + solution_values.push_back (dof->second); + } + solution.set (indices, solution_values); + + // now also set appropriate + // values for the rhs + for (unsigned int i=0; i<solution_values.size(); ++i) + solution_values[i] *= matrix.diag_element(indices[i]); + + right_hand_side.set (indices, solution_values); + } + else + { + // clear_rows() is a collective operation so we still have to call + // it: + std::vector<types::global_dof_index> constrained_rows; + matrix.clear_rows (constrained_rows, 1.); + } + + // clean up + matrix.compress (VectorOperation::insert); + solution.compress (VectorOperation::insert); + right_hand_side.compress (VectorOperation::insert); + } + + + + template <typename TrilinosMatrix, typename TrilinosBlockVector> + void + apply_block_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + TrilinosMatrix &matrix, + TrilinosBlockVector &solution, + TrilinosBlockVector &right_hand_side, + const bool eliminate_columns) + { + Assert (eliminate_columns == false, ExcNotImplemented()); + + Assert (matrix.n() == right_hand_side.size(), + ExcDimensionMismatch(matrix.n(), right_hand_side.size())); + Assert (matrix.n() == solution.size(), + ExcDimensionMismatch(matrix.n(), solution.size())); + Assert (matrix.n_block_rows() == matrix.n_block_cols(), + ExcNotQuadratic()); + + const unsigned int n_blocks = matrix.n_block_rows(); + + // We need to find the subdivision + // into blocks for the boundary values. + // To this end, generate a vector of + // maps with the respective indices. + std::vector<std::map<types::global_dof_index,double> > block_boundary_values(n_blocks); + { + int block=0; + types::global_dof_index offset = 0; + for (std::map<types::global_dof_index,double>::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + { + if (dof->first >= matrix.block(block,0).m() + offset) + { + offset += matrix.block(block,0).m(); + block++; + } + const types::global_dof_index index = dof->first - offset; + block_boundary_values[block].insert( + std::pair<types::global_dof_index, double> (index,dof->second)); + } + } + + // Now call the non-block variants on + // the diagonal subblocks and the + // solution/rhs. + for (unsigned int block=0; block<n_blocks; ++block) + TrilinosWrappers::apply_boundary_values(block_boundary_values[block], + matrix.block(block,block), + solution.block(block), + right_hand_side.block(block), + eliminate_columns); + + // Finally, we need to do something + // about the off-diagonal matrices. This + // is luckily not difficult. Just clear + // the whole row. + for (unsigned int block_m=0; block_m<n_blocks; ++block_m) + { + const std::pair<types::global_dof_index, types::global_dof_index> local_range + = matrix.block(block_m,0).local_range(); + + std::vector<types::global_dof_index> constrained_rows; + for (std::map<types::global_dof_index,double>::const_iterator + dof = block_boundary_values[block_m].begin(); + dof != block_boundary_values[block_m].end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + constrained_rows.push_back (dof->first); + + for (unsigned int block_n=0; block_n<n_blocks; ++block_n) + if (block_m != block_n) + matrix.block(block_m,block_n).clear_rows(constrained_rows); + } + } + } +} + + + + +void +apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + TrilinosWrappers::SparseMatrix &matrix, + TrilinosWrappers::Vector &solution, + TrilinosWrappers::Vector &right_hand_side, + const bool eliminate_columns) +{ + // simply redirect to the generic function + // used for both trilinos matrix types + internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution, + right_hand_side, eliminate_columns); +} + + + +void +apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + TrilinosWrappers::SparseMatrix &matrix, + TrilinosWrappers::MPI::Vector &solution, + TrilinosWrappers::MPI::Vector &right_hand_side, + const bool eliminate_columns) +{ + // simply redirect to the generic function + // used for both trilinos matrix types + internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution, + right_hand_side, eliminate_columns); +} + + + +void +apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + TrilinosWrappers::BlockSparseMatrix &matrix, + TrilinosWrappers::BlockVector &solution, + TrilinosWrappers::BlockVector &right_hand_side, + const bool eliminate_columns) +{ + internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix, + solution, right_hand_side, + eliminate_columns); +} + + + +void +apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values, + TrilinosWrappers::BlockSparseMatrix &matrix, + TrilinosWrappers::MPI::BlockVector &solution, + TrilinosWrappers::MPI::BlockVector &right_hand_side, + const bool eliminate_columns) +{ + internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix, + solution, right_hand_side, + eliminate_columns); +} + +#endif + + +DEAL_II_NAMESPACE_CLOSE -- 2.39.5