From fad0dfb7a9412258ee4a2f8910bebbf6c55316ed Mon Sep 17 00:00:00 2001 From: Matthias Maier 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_WITH_PETSC +# include +# include +# include +# include +# include +#endif + +#ifdef DEAL_II_WITH_TRILINOS +# include +# include +# include +# include +#endif + +#include + + +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace MatrixCreator +{ + namespace internal + { + /** + * Convenience abbreviation for + * pairs of DoF handler cell + * iterators. This type works + * just like a + * std::pair + * but is templatized on the + * dof handler that should be used. + */ + template + 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 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 + inline + IteratorRange:: + IteratorRange (const active_cell_iterator &first, + const active_cell_iterator &second) + : + first (first), + second (second) + {} + + + + template + inline + IteratorRange::IteratorRange (const iterator_pair &ip) + : + first (ip.first), + second (ip.second) + {} + + + + namespace AssemblerData + { + template + struct Scratch + { + Scratch (const ::dealii::hp::FECollection &fe, + const UpdateFlags update_flags, + const Function *coefficient, + const Function *rhs_function, + const ::dealii::hp::QCollection &quadrature, + const ::dealii::hp::MappingCollection &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 (fe_collection.n_components())), + rhs_values(quadrature_collection.max_n_quadrature_points()), + rhs_vector_values(quadrature_collection.max_n_quadrature_points(), + dealii::Vector (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 &fe_collection; + const ::dealii::hp::QCollection &quadrature_collection; + const ::dealii::hp::MappingCollection &mapping_collection; + + ::dealii::hp::FEValues x_fe_values; + + std::vector coefficient_values; + std::vector > coefficient_vector_values; + std::vector rhs_values; + std::vector > rhs_vector_values; + + const Function *coefficient; + const Function *rhs_function; + + const UpdateFlags update_flags; + }; + + + template + struct CopyData + { + std::vector dof_indices; + FullMatrix cell_matrix; + dealii::Vector cell_rhs; + const ConstraintMatrix *constraints; + }; + } + + + template + void mass_assembler (const CellIterator &cell, + MatrixCreator::internal::AssemblerData::Scratch &data, + MatrixCreator::internal::AssemblerData::CopyData ©_data) + { + data.x_fe_values.reinit (cell); + const FEValues &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 &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(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(n_components)); + data.coefficient->vector_value_list (fe_values.get_quadrature_points(), + data.coefficient_vector_values); + } + } + + + double add_data; + const std::vector &JxW = fe_values.get_JxW_values(); + for (unsigned int i=0; in_components==1) + for (unsigned int point=0; pointn_components==1) + for (unsigned int point=0; pointn_components==1) + for (unsigned int point=0; pointn_components==1) + for (unsigned int point=0; point + void laplace_assembler (const CellIterator &cell, + MatrixCreator::internal::AssemblerData::Scratch &data, + MatrixCreator::internal::AssemblerData::CopyData ©_data) + { + data.x_fe_values.reinit (cell); + const FEValues &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 &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(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(n_components)); + data.coefficient->vector_value_list (fe_values.get_quadrature_points(), + data.coefficient_vector_values); + } + } + + + const std::vector &JxW = fe_values.get_JxW_values(); + double add_data; + for (unsigned int i=0; i *grad_phi_i = + &fe_values.shape_grad(i,0); + + // can use symmetry + for (unsigned int j=i; j *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; pointn_components==1) + for (unsigned int point=0; pointn_components==1) + for (unsigned int point=0; pointn_components==1) + for (unsigned int point=0; point + void copy_local_to_global (const AssemblerData::CopyData &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 + struct CopyData + { + CopyData() {}; + + CopyData(CopyData const &data); + + unsigned int dofs_per_cell; + std::vector dofs; + std::vector > dof_is_on_face; + typename DoFHandlerType::active_cell_iterator cell; + std::vector > cell_matrix; + std::vector > cell_vector; + }; + + template + CopyData::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 + void create_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function *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 fe_collection (dof.get_fe()); + hp::QCollection q_collection (q); + hp::MappingCollection mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch + 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 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::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler::active_cell_iterator,number>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx11::_1, &matrix, (Vector *)0), + assembler_data, + copy_data); + } + + + + template + void create_mass_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function *const coefficient, + const ConstraintMatrix &constraints) + { + create_mass_matrix(StaticMappingQ1::mapping, dof, + q, matrix, coefficient, constraints); + } + + + + template + void create_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function *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 fe_collection (dof.get_fe()); + hp::QCollection q_collection (q); + hp::MappingCollection mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (fe_collection, + update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q_collection, mapping_collection); + MatrixCreator::internal::AssemblerData::CopyData 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::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler::active_cell_iterator,number>, + std_cxx11::bind(&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx11::_1, &matrix, &rhs_vector), + assembler_data, + copy_data); + } + + + + template + void create_mass_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function *const coefficient, + const ConstraintMatrix &constraints) + { + create_mass_matrix(StaticMappingQ1::mapping, + dof, q, matrix, rhs, rhs_vector, coefficient, + constraints); + } + + + + template + void create_mass_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function *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 + 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 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::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler::active_cell_iterator,number>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx11::_1, &matrix, (Vector *)0), + assembler_data, + copy_data); + } + + + + template + void create_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function *const coefficient, + const ConstraintMatrix &constraints) + { + create_mass_matrix(hp::StaticMappingQ1::mapping_collection, + dof, q, matrix, coefficient, constraints); + } + + + + template + void create_mass_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function *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 + assembler_data (dof.get_fe(), + update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData 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::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler::active_cell_iterator,number>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx11::_1, &matrix, &rhs_vector), + assembler_data, + copy_data); + } + + + + template + void create_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function *const coefficient, + const ConstraintMatrix &constraints) + { + create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, + matrix, rhs, rhs_vector, coefficient, constraints); + } + + + + namespace internal + { + template + void + static inline + create_boundary_mass_matrix_1 (typename DoFHandler::active_cell_iterator const &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData > ©_data, + Mapping const &mapping, + FiniteElement const &fe, + Quadrature const &q, + typename FunctionMap::type const &boundary_functions, + Function const *const coefficient, + std::vector 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 fe_values (mapping, fe, q, update_flags); + + // two variables for the coefficient, + // one for the two cases indicated in + // the name + std::vector coefficient_values (fe_values.n_quadrature_points, 1.); + std::vector > coefficient_vector_values (fe_values.n_quadrature_points, + Vector(n_components)); + const bool coefficient_is_vector = (coefficient != 0 && coefficient->n_components != 1) + ? true : false; + + std::vector rhs_values_scalar (fe_values.n_quadrature_points); + std::vector > rhs_values_system (fe_values.n_quadrature_points, + Vector(n_function_components)); + + copy_data.dofs.resize(copy_data.dofs_per_cell); + cell->get_dof_indices (copy_data.dofs); + + std::vector 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::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 (copy_data.dofs_per_cell, + copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector (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 > normal_adjustment(fe_values.n_quadrature_points, + std::vector(n_components, 1.)); + + for (unsigned int comp = 0; comp &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::H1) && + base.conforms(FiniteElementData::Hdiv)) + for (unsigned int point=0; pointface(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; pointface(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 (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 + void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData > const ©_data, + typename FunctionMap::type const &boundary_functions, + std::vector const &dof_to_boundary_mapping, + SparseMatrix &matrix, + Vector &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::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 + 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 > &/*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 const &/*component_mapping*/) + { + Assert(false,ExcNotImplemented()); + } + } + + + + + template + void + create_boundary_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const typename FunctionMap::type &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function *const coefficient, + std::vector component_mapping) + { + // what would that be in 1d? the + // identity matrix on the boundary + // dofs? + if (dim == 1) + { + Assert (false, ExcNotImplemented()); + return; + } + + const FiniteElement &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 > copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > + (std_cxx11::bind(&internal::create_boundary_mass_matrix_1,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 > const &)> > (std_cxx11::bind( + &internal::copy_boundary_mass_matrix_1, + 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 + void + create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler::active_cell_iterator const + &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary + ::CopyData > ©_data, + hp::MappingCollection const &mapping, + hp::FECollection const &fe_collection, + hp::QCollection const &q, + const typename FunctionMap::type &boundary_functions, + Function const *const coefficient, + std::vector 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 &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 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 coefficient_values; + std::vector > coefficient_vector_values; + + std::vector rhs_values_scalar; + std::vector > rhs_values_system; + + std::vector 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::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 &fe_values = x_fe_values.get_present_fe_values (); + + copy_data.cell_matrix.push_back(FullMatrix (copy_data.dofs_per_cell, + copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); + + if (fe_is_system) + // FE has several components + { + rhs_values_system.resize (fe_values.n_quadrature_points, + Vector(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(n_components)); + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_vector_values); + for (unsigned int point=0; pointface(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; pointface(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 (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 + void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary + ::CopyData > const ©_data, + typename FunctionMap::type const &boundary_functions, + std::vector const &dof_to_boundary_mapping, + SparseMatrix &matrix, + Vector &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::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(0); + for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != hp::DoFHandler::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 max_diag_entry) + max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i)); +#endif + + for (unsigned int i=0; i + void create_boundary_mass_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const typename FunctionMap::type &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function *const a, + std::vector component_mapping) + { + create_boundary_mass_matrix(StaticMappingQ1::mapping, dof, q, + matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); + } + + + + template + void + create_boundary_mass_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const typename FunctionMap::type &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function *const coefficient, + std::vector 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 &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 > copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > + (std_cxx11::bind( &create_hp_boundary_mass_matrix_1,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 > const &)> > ( + std_cxx11::bind( ©_hp_boundary_mass_matrix_1, + 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 + void create_boundary_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const typename FunctionMap::type &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function *const a, + std::vector component_mapping) + { + create_boundary_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, + matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); + } + + + + template + void create_laplace_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function *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 fe_collection (dof.get_fe()); + hp::QCollection q_collection (q); + hp::MappingCollection mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch + 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 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::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx11::_1, + &matrix, + (Vector *)NULL), + assembler_data, + copy_data); + } + + + + template + void create_laplace_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function *const coefficient, + const ConstraintMatrix &constraints) + { + create_laplace_matrix(StaticMappingQ1::mapping, + dof, q, matrix, coefficient, constraints); + } + + + + template + void create_laplace_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function *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 fe_collection (dof.get_fe()); + hp::QCollection q_collection (q); + hp::MappingCollection mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (fe_collection, + update_gradients | update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q_collection, mapping_collection); + MatrixCreator::internal::AssemblerData::CopyData 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::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx11::_1, + &matrix, + &rhs_vector), + assembler_data, + copy_data); + } + + + + template + void create_laplace_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function *const coefficient, + const ConstraintMatrix &constraints) + { + create_laplace_matrix(StaticMappingQ1::mapping, dof, q, + matrix, rhs, rhs_vector, coefficient, constraints); + } + + + + template + void create_laplace_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function *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 + 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 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::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx11::_1, + &matrix, + (Vector *)0), + assembler_data, + copy_data); + } + + + + template + void create_laplace_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function *const coefficient, + const ConstraintMatrix &constraints) + { + create_laplace_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, + matrix, coefficient, constraints); + } + + + + template + void create_laplace_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function *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 + assembler_data (dof.get_fe(), + update_gradients | update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData 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::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, + std_cxx11::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx11::_1, + &matrix, + &rhs_vector), + assembler_data, + copy_data); + } + + + + template + void create_laplace_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function *const coefficient, + const ConstraintMatrix &constraints) + { + create_laplace_matrix(hp::StaticMappingQ1::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_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 + (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + + template + void MatrixCreator::create_boundary_mass_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const FunctionMap::type &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const a, + std::vector); + + template + void MatrixCreator::create_boundary_mass_matrix + (const Mapping &, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const FunctionMap::type &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const a, + std::vector); + + template + void + MatrixCreator::create_boundary_mass_matrix + (const hp::MappingCollection&, + const hp::DoFHandler&, + const hp::QCollection&, + SparseMatrix&, + const FunctionMap::type&, + Vector&, + std::vector&, + const Function * const, + std::vector); + + template + void MatrixCreator::create_boundary_mass_matrix + (const hp::DoFHandler&, + const hp::QCollection&, + SparseMatrix&, + const FunctionMap::type&, + Vector&, + std::vector&, + const Function * const, + std::vector); + +// same for float + template + void MatrixCreator::create_mass_matrix + (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_mass_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + +#endif + } + + +//TODO[SP]: replace by +// 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 + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + +#endif + +#if deal_II_dimension == deal_II_space_dimension + + + template + void MatrixCreator::create_mass_matrix + (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix + (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + +// non-hp versions of create_laplace_matrix + template + void MatrixCreator::create_laplace_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix + (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix + (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix + (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + +// hp versions of create_laplace_matrix + template + void MatrixCreator::create_laplace_matrix + (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + template + void MatrixCreator::create_laplace_matrix + (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + +#endif + +// same for float +#if deal_II_dimension <= deal_II_space_dimension + + template + void MatrixCreator::create_mass_matrix + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + +#endif + +#if deal_II_dimension == deal_II_space_dimension + + + template + void MatrixCreator::create_mass_matrix + (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient, + const ConstraintMatrix &constraints); + + template + void MatrixCreator::create_mass_matrix + (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * 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 - * std::pair - * but is templatized on the - * dof handler that should be used. - */ - template - 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 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 - inline - IteratorRange:: - IteratorRange (const active_cell_iterator &first, - const active_cell_iterator &second) - : - first (first), - second (second) - {} - - - - template - inline - IteratorRange::IteratorRange (const iterator_pair &ip) - : - first (ip.first), - second (ip.second) - {} - - - - namespace AssemblerData - { - template - struct Scratch - { - Scratch (const ::dealii::hp::FECollection &fe, - const UpdateFlags update_flags, - const Function *coefficient, - const Function *rhs_function, - const ::dealii::hp::QCollection &quadrature, - const ::dealii::hp::MappingCollection &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 (fe_collection.n_components())), - rhs_values(quadrature_collection.max_n_quadrature_points()), - rhs_vector_values(quadrature_collection.max_n_quadrature_points(), - dealii::Vector (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 &fe_collection; - const ::dealii::hp::QCollection &quadrature_collection; - const ::dealii::hp::MappingCollection &mapping_collection; - - ::dealii::hp::FEValues x_fe_values; - - std::vector coefficient_values; - std::vector > coefficient_vector_values; - std::vector rhs_values; - std::vector > rhs_vector_values; - - const Function *coefficient; - const Function *rhs_function; - - const UpdateFlags update_flags; - }; - - - template - struct CopyData - { - std::vector dof_indices; - FullMatrix cell_matrix; - dealii::Vector cell_rhs; - const ConstraintMatrix *constraints; - }; - } - - - template - void mass_assembler (const CellIterator &cell, - MatrixCreator::internal::AssemblerData::Scratch &data, - MatrixCreator::internal::AssemblerData::CopyData ©_data) - { - data.x_fe_values.reinit (cell); - const FEValues &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 &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(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(n_components)); - data.coefficient->vector_value_list (fe_values.get_quadrature_points(), - data.coefficient_vector_values); - } - } - - - double add_data; - const std::vector &JxW = fe_values.get_JxW_values(); - for (unsigned int i=0; in_components==1) - for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; point - void laplace_assembler (const CellIterator &cell, - MatrixCreator::internal::AssemblerData::Scratch &data, - MatrixCreator::internal::AssemblerData::CopyData ©_data) - { - data.x_fe_values.reinit (cell); - const FEValues &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 &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(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(n_components)); - data.coefficient->vector_value_list (fe_values.get_quadrature_points(), - data.coefficient_vector_values); - } - } - - - const std::vector &JxW = fe_values.get_JxW_values(); - double add_data; - for (unsigned int i=0; i *grad_phi_i = - &fe_values.shape_grad(i,0); - - // can use symmetry - for (unsigned int j=i; j *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; pointn_components==1) - for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; pointn_components==1) - for (unsigned int point=0; point - void copy_local_to_global (const AssemblerData::CopyData &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 - struct CopyData - { - CopyData() {}; - - CopyData(CopyData const &data); - - unsigned int dofs_per_cell; - std::vector dofs; - std::vector > dof_is_on_face; - typename DoFHandlerType::active_cell_iterator cell; - std::vector > cell_matrix; - std::vector > cell_vector; - }; - - template - CopyData::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 - void create_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function *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 fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - hp::MappingCollection mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch - 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 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::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::mass_assembler::active_cell_iterator,number>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global, Vector >, - std_cxx11::_1, &matrix, (Vector *)0), - assembler_data, - copy_data); - } - - - - template - void create_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function *const coefficient, - const ConstraintMatrix &constraints) - { - create_mass_matrix(StaticMappingQ1::mapping, dof, - q, matrix, coefficient, constraints); - } - - - - template - void create_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function *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 fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - hp::MappingCollection mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch - assembler_data (fe_collection, - update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q_collection, mapping_collection); - MatrixCreator::internal::AssemblerData::CopyData 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::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::mass_assembler::active_cell_iterator,number>, - std_cxx11::bind(&MatrixCreator::internal:: - copy_local_to_global, Vector >, - std_cxx11::_1, &matrix, &rhs_vector), - assembler_data, - copy_data); - } - - - - template - void create_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function *const coefficient, - const ConstraintMatrix &constraints) - { - create_mass_matrix(StaticMappingQ1::mapping, - dof, q, matrix, rhs, rhs_vector, coefficient, - constraints); - } - - - - template - void create_mass_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function *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 - 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 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::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::mass_assembler::active_cell_iterator,number>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global, Vector >, - std_cxx11::_1, &matrix, (Vector *)0), - assembler_data, - copy_data); - } - - - - template - void create_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function *const coefficient, - const ConstraintMatrix &constraints) - { - create_mass_matrix(hp::StaticMappingQ1::mapping_collection, - dof, q, matrix, coefficient, constraints); - } - - - - template - void create_mass_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function *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 - assembler_data (dof.get_fe(), - update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q, mapping); - MatrixCreator::internal::AssemblerData::CopyData 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::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::mass_assembler::active_cell_iterator,number>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global, Vector >, - std_cxx11::_1, &matrix, &rhs_vector), - assembler_data, - copy_data); - } - - - - template - void create_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function *const coefficient, - const ConstraintMatrix &constraints) - { - create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, - matrix, rhs, rhs_vector, coefficient, constraints); - } - - - - namespace - { - template - void - create_boundary_mass_matrix_1 (typename DoFHandler::active_cell_iterator const &cell, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData > ©_data, - Mapping const &mapping, - FiniteElement const &fe, - Quadrature const &q, - typename FunctionMap::type const &boundary_functions, - Function const *const coefficient, - std::vector 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 fe_values (mapping, fe, q, update_flags); - - // two variables for the coefficient, - // one for the two cases indicated in - // the name - std::vector coefficient_values (fe_values.n_quadrature_points, 1.); - std::vector > coefficient_vector_values (fe_values.n_quadrature_points, - Vector(n_components)); - const bool coefficient_is_vector = (coefficient != 0 && coefficient->n_components != 1) - ? true : false; - - std::vector rhs_values_scalar (fe_values.n_quadrature_points); - std::vector > rhs_values_system (fe_values.n_quadrature_points, - Vector(n_function_components)); - - copy_data.dofs.resize(copy_data.dofs_per_cell); - cell->get_dof_indices (copy_data.dofs); - - std::vector 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::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 (copy_data.dofs_per_cell, - copy_data.dofs_per_cell)); - copy_data.cell_vector.push_back(Vector (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 > normal_adjustment(fe_values.n_quadrature_points, - std::vector(n_components, 1.)); - - for (unsigned int comp = 0; comp &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::H1) && - base.conforms(FiniteElementData::Hdiv)) - for (unsigned int point=0; pointface(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; pointface(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 (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 - void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData > const ©_data, - typename FunctionMap::type const &boundary_functions, - std::vector const &dof_to_boundary_mapping, - SparseMatrix &matrix, - Vector &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::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 - 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 > &/*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 const &/*component_mapping*/) - { - Assert(false,ExcNotImplemented()); - } - } - - - - - template - void - create_boundary_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const coefficient, - std::vector component_mapping) - { - // what would that be in 1d? the - // identity matrix on the boundary - // dofs? - if (dim == 1) - { - Assert (false, ExcNotImplemented()); - return; - } - - const FiniteElement &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 > copy_data; - - WorkStream::run(dof.begin_active(),dof.end(), - static_cast::active_cell_iterator - const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > - (std_cxx11::bind(&create_boundary_mass_matrix_1,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 > const &)> > (std_cxx11::bind( - ©_boundary_mass_matrix_1, - 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 - void - create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler::active_cell_iterator const - &cell, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary - ::CopyData > ©_data, - hp::MappingCollection const &mapping, - hp::FECollection const &fe_collection, - hp::QCollection const &q, - const typename FunctionMap::type &boundary_functions, - Function const *const coefficient, - std::vector 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 &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 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 coefficient_values; - std::vector > coefficient_vector_values; - - std::vector rhs_values_scalar; - std::vector > rhs_values_system; - - std::vector 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::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 &fe_values = x_fe_values.get_present_fe_values (); - - copy_data.cell_matrix.push_back(FullMatrix (copy_data.dofs_per_cell, - copy_data.dofs_per_cell)); - copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); - - if (fe_is_system) - // FE has several components - { - rhs_values_system.resize (fe_values.n_quadrature_points, - Vector(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(n_components)); - coefficient->vector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - for (unsigned int point=0; pointface(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; pointface(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 (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 - void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary - ::CopyData > const ©_data, - typename FunctionMap::type const &boundary_functions, - std::vector const &dof_to_boundary_mapping, - SparseMatrix &matrix, - Vector &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::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(0); - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != hp::DoFHandler::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 max_diag_entry) - max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i)); -#endif - - for (unsigned int i=0; i - void create_boundary_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const typename FunctionMap::type &rhs, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const a, - std::vector component_mapping) - { - create_boundary_mass_matrix(StaticMappingQ1::mapping, dof, q, - matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); - } - - - - template - void - create_boundary_mass_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const coefficient, - std::vector 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 &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 > copy_data; - - WorkStream::run(dof.begin_active(),dof.end(), - static_cast::active_cell_iterator - const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > - (std_cxx11::bind( &create_hp_boundary_mass_matrix_1,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 > const &)> > ( - std_cxx11::bind( ©_hp_boundary_mass_matrix_1, - 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 - void create_boundary_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const typename FunctionMap::type &rhs, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const a, - std::vector component_mapping) - { - create_boundary_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, - matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); - } - - - - template - void create_laplace_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function *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 fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - hp::MappingCollection mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch - 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 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::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global, Vector >, - std_cxx11::_1, - &matrix, - (Vector *)NULL), - assembler_data, - copy_data); - } - - - - template - void create_laplace_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function *const coefficient, - const ConstraintMatrix &constraints) - { - create_laplace_matrix(StaticMappingQ1::mapping, - dof, q, matrix, coefficient, constraints); - } - - - - template - void create_laplace_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function *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 fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - hp::MappingCollection mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch - assembler_data (fe_collection, - update_gradients | update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q_collection, mapping_collection); - MatrixCreator::internal::AssemblerData::CopyData 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::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global, Vector >, - std_cxx11::_1, - &matrix, - &rhs_vector), - assembler_data, - copy_data); - } - - - - template - void create_laplace_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function *const coefficient, - const ConstraintMatrix &constraints) - { - create_laplace_matrix(StaticMappingQ1::mapping, dof, q, - matrix, rhs, rhs_vector, coefficient, constraints); - } - - - - template - void create_laplace_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function *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 - 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 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::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global, Vector >, - std_cxx11::_1, - &matrix, - (Vector *)0), - assembler_data, - copy_data); - } - - - - template - void create_laplace_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function *const coefficient, - const ConstraintMatrix &constraints) - { - create_laplace_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, - matrix, coefficient, constraints); - } - - - - template - void create_laplace_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function *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 - assembler_data (dof.get_fe(), - update_gradients | update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q, mapping); - MatrixCreator::internal::AssemblerData::CopyData 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::active_cell_iterator>(dof.end()), - &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, - std_cxx11::bind (&MatrixCreator::internal:: - copy_local_to_global, Vector >, - std_cxx11::_1, - &matrix, - &rhs_vector), - assembler_data, - copy_data); - } - - - - template - void create_laplace_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function *const coefficient, - const ConstraintMatrix &constraints) - { - create_laplace_matrix(hp::StaticMappingQ1::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 - void - apply_boundary_values (const std::map &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 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 constrained_rows; - for (std::map::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 indices; - std::vector solution_values; - for (std::map::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 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 &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 &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 &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 > block_boundary_values(n_blocks); - { - int block = 0; - dealii::types::global_dof_index offset = 0; - for (std::map::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 (index,dof->second)); - } - } - - // Now call the non-block variants on - // the diagonal subblocks and the - // solution/rhs. - for (unsigned int block=0; block local_range - = matrix.block(block_m,0).local_range(); - - std::vector constrained_rows; - for (std::map::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 - void - apply_boundary_values (const std::map &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 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 constrained_rows; - for (std::map::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 indices; - std::vector solution_values; - for (std::map::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 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 - void - apply_block_boundary_values (const std::map &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 > block_boundary_values(n_blocks); - { - int block=0; - types::global_dof_index offset = 0; - for (std::map::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 (index,dof->second)); - } - } - - // Now call the non-block variants on - // the diagonal subblocks and the - // solution/rhs. - for (unsigned int block=0; block local_range - = matrix.block(block_m,0).local_range(); - - std::vector constrained_rows; - for (std::map::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 &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 &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 &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 &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 (const std::map &boundary_values, - SparseMatrix &matrix, - Vector &solution, - Vector &right_hand_side, - const bool eliminate_columns); - template - void - apply_boundary_values (const std::map &boundary_values, - SparseMatrix &matrix, - Vector &solution, - Vector &right_hand_side, - const bool eliminate_columns); - - template - void - apply_boundary_values (const std::map &boundary_values, - BlockSparseMatrix &matrix, - BlockVector &solution, - BlockVector &right_hand_side, - const bool eliminate_columns); - template - void - apply_boundary_values (const std::map &boundary_values, - BlockSparseMatrix &matrix, - BlockVector &solution, - BlockVector &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 - (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - - - template - void MatrixCreator::create_boundary_mass_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const FunctionMap::type &rhs, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const a, - std::vector); - - template - void MatrixCreator::create_boundary_mass_matrix - (const Mapping &, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const FunctionMap::type &rhs, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const a, - std::vector); - - template - void - MatrixCreator::create_boundary_mass_matrix - (const hp::MappingCollection&, - const hp::DoFHandler&, - const hp::QCollection&, - SparseMatrix&, - const FunctionMap::type&, - Vector&, - std::vector&, - const Function * const, - std::vector); - - template - void MatrixCreator::create_boundary_mass_matrix - (const hp::DoFHandler&, - const hp::QCollection&, - SparseMatrix&, - const FunctionMap::type&, - Vector&, - std::vector&, - const Function * const, - std::vector); - -// same for float - template - void MatrixCreator::create_mass_matrix - (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - -#endif - } - -//TODO[SP]: replace by -// 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 - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - - -#endif - -#if deal_II_dimension == deal_II_space_dimension + void MatrixTools::apply_boundary_values + (const std::map &boundary_values, + SparseMatrix &matrix, + Vector &solution, + Vector &right_hand_side, + const bool eliminate_columns); - template - void MatrixCreator::create_mass_matrix - (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - - template - void MatrixCreator::create_mass_matrix - (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - -// non-hp versions of create_laplace_matrix - template - void MatrixCreator::create_laplace_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix - (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix - (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - -// hp versions of create_laplace_matrix - template - void MatrixCreator::create_laplace_matrix - (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_laplace_matrix - (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - -#endif - -// same for float -#if deal_II_dimension <= deal_II_space_dimension - - template - void MatrixCreator::create_mass_matrix - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - - template - void MatrixCreator::create_mass_matrix - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - - -#endif - -#if deal_II_dimension == deal_II_space_dimension - - - template - void MatrixCreator::create_mass_matrix - (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - - template - void MatrixCreator::create_mass_matrix - (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - -#endif + void MatrixTools::apply_boundary_values + (const std::map &boundary_values, + BlockSparseMatrix &matrix, + BlockVector &solution, + BlockVector &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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_WITH_PETSC +# include +# include +# include +# include +# include +#endif + +#ifdef DEAL_II_WITH_TRILINOS +# include +# include +# include +# include +#endif + +#include + + +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +#ifdef DEAL_II_WITH_PETSC + +namespace internal +{ + namespace PETScWrappers + { + template + void + apply_boundary_values (const std::map &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 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 constrained_rows; + for (std::map::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 indices; + std::vector solution_values; + for (std::map::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 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 &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 &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 &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 > block_boundary_values(n_blocks); + { + int block = 0; + dealii::types::global_dof_index offset = 0; + for (std::map::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 (index,dof->second)); + } + } + + // Now call the non-block variants on + // the diagonal subblocks and the + // solution/rhs. + for (unsigned int block=0; block local_range + = matrix.block(block_m,0).local_range(); + + std::vector constrained_rows; + for (std::map::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 + void + apply_boundary_values (const std::map &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 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 constrained_rows; + for (std::map::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 indices; + std::vector solution_values; + for (std::map::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 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 + void + apply_block_boundary_values (const std::map &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 > block_boundary_values(n_blocks); + { + int block=0; + types::global_dof_index offset = 0; + for (std::map::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 (index,dof->second)); + } + } + + // Now call the non-block variants on + // the diagonal subblocks and the + // solution/rhs. + for (unsigned int block=0; block local_range + = matrix.block(block_m,0).local_range(); + + std::vector constrained_rows; + for (std::map::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 &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 &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 &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 &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