From 1e171e3be4cf80a9d01b1345024092ca8758a289 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Mon, 29 Feb 2016 11:22:11 -0500 Subject: [PATCH] move to a .templates.h this avoid warnings with clang about undefined functions --- .../numerics/matrix_creator.templates.h | 1937 +++++++++++++++++ source/numerics/matrix_creator.cc | 1908 +--------------- 2 files changed, 1938 insertions(+), 1907 deletions(-) create mode 100644 include/deal.II/numerics/matrix_creator.templates.h 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..4541403de0 --- /dev/null +++ b/include/deal.II/numerics/matrix_creator.templates.h @@ -0,0 +1,1937 @@ +// --------------------------------------------------------------------- +// +// 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); + } + } + + + number 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,number > ©_data, + Mapping const &mapping, + FiniteElement const &fe, + Quadrature const &q, + std::map*> 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,number > const ©_data, + std::map*> 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,float> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData,float > &/*copy_data*/, + Mapping<1,3> const &, + FiniteElement<1,3> const &, + Quadrature<0> const &, + std::map*> const &/*boundary_functions*/, + Function<3,float> const *const /*coefficient*/, + std::vector const &/*component_mapping*/) + { + Assert(false,ExcNotImplemented()); + } + + template <> + void + inline + create_boundary_mass_matrix_1<1,3,double> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData,double > &/*copy_data*/, + Mapping<1,3> const &, + FiniteElement<1,3> const &, + Quadrature<0> const &, + std::map*> const &/*boundary_functions*/, + Function<3,double> const *const /*coefficient*/, + std::vector const &/*component_mapping*/) + { + Assert(false,ExcNotImplemented()); + } + + template <> + void + inline + create_boundary_mass_matrix_1<1,3,long double> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData,long double > &/*copy_data*/, + Mapping<1,3> const &, + FiniteElement<1,3> const &, + Quadrature<0> const &, + std::map*> const &/*boundary_functions*/, + Function<3,long double> 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 std::map*> &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,number > copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData,number > &)> > + (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,number> 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,number > ©_data, + hp::MappingCollection const &mapping, + hp::FECollection const &fe_collection, + hp::QCollection const &q, + const std::map*> &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,number > const ©_data, + std::map*> 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::abs(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 std::map*> &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 std::map*> &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,number > copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData,number > &)> > + (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,number > 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 std::map*> &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/matrix_creator.cc b/source/numerics/matrix_creator.cc index 6f581d9c99..1bb6205054 100644 --- a/source/numerics/matrix_creator.cc +++ b/source/numerics/matrix_creator.cc @@ -13,1915 +13,10 @@ // // --------------------------------------------------------------------- -#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 +#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); - } - } - - - number 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,number > ©_data, - Mapping const &mapping, - FiniteElement const &fe, - Quadrature const &q, - std::map*> 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,number > const ©_data, - std::map*> 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,float> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData,float > &/*copy_data*/, - Mapping<1,3> const &, - FiniteElement<1,3> const &, - Quadrature<0> const &, - std::map*> const &/*boundary_functions*/, - Function<3,float> const *const /*coefficient*/, - std::vector const &/*component_mapping*/) - { - Assert(false,ExcNotImplemented()); - } - - template <> - void - create_boundary_mass_matrix_1<1,3,double> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData,double > &/*copy_data*/, - Mapping<1,3> const &, - FiniteElement<1,3> const &, - Quadrature<0> const &, - std::map*> const &/*boundary_functions*/, - Function<3,double> const *const /*coefficient*/, - std::vector const &/*component_mapping*/) - { - Assert(false,ExcNotImplemented()); - } - - template <> - void - create_boundary_mass_matrix_1<1,3,long double> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData,long double > &/*copy_data*/, - Mapping<1,3> const &, - FiniteElement<1,3> const &, - Quadrature<0> const &, - std::map*> const &/*boundary_functions*/, - Function<3,long double> 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 std::map*> &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,number > copy_data; - - WorkStream::run(dof.begin_active(),dof.end(), - static_cast::active_cell_iterator - const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData,number > &)> > - (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,number> 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,number > ©_data, - hp::MappingCollection const &mapping, - hp::FECollection const &fe_collection, - hp::QCollection const &q, - const std::map*> &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,number > const ©_data, - std::map*> 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::abs(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 std::map*> &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 std::map*> &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,number > copy_data; - - WorkStream::run(dof.begin_active(),dof.end(), - static_cast::active_cell_iterator - const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData,number > &)> > - (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,number > 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 std::map*> &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 - // explicit instantiations #define SPLIT_INSTANTIATIONS_COUNT 3 @@ -1930,5 +25,4 @@ namespace MatrixCreator #endif #include "matrix_creator.inst" - DEAL_II_NAMESPACE_CLOSE -- 2.39.5