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