#include <base/quadrature.h>
#include <fe/fe.h>
#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
+#include <fe/mapping_collection.h>
#include <numerics/matrices.h>
#include <lac/vector.h>
#include <lac/block_vector.h>
DEAL_II_NAMESPACE_OPEN
-template <int dim>
+template <typename DH>
inline
-MatrixCreator::IteratorRange<dim>::
+MatrixCreator::IteratorRange<DH>::
IteratorRange (const active_cell_iterator &first,
- const active_cell_iterator &second)
+ const active_cell_iterator &second)
:
first (first),
second (second)
-template <int dim>
+template <typename DH>
inline
-MatrixCreator::IteratorRange<dim>::IteratorRange (const iterator_pair &ip)
+MatrixCreator::IteratorRange<DH>::IteratorRange (const iterator_pair &ip)
:
first (ip.first),
second (ip.second)
const Quadrature<dim> &q,
SparseMatrix<number> &matrix,
const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
+ const IteratorRange<DoFHandler<dim> > range,
Threads::ThreadMutex &mutex);
create_mass_matrix_1_t p = &MatrixCreator::template create_mass_matrix_1<dim,number>;
for (unsigned int thread=0; thread<n_threads; ++thread)
const Quadrature<dim> &q,
SparseMatrix<number> &matrix,
const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
+ const IteratorRange<DoFHandler<dim> > range,
Threads::ThreadMutex &mutex)
{
UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values);
const Function<dim> &rhs,
Vector<double> &rhs_vector,
const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
+ const IteratorRange<DoFHandler<dim> > range,
Threads::ThreadMutex &mutex);
create_mass_matrix_2_t p = &MatrixCreator::template create_mass_matrix_2<dim,number>;
for (unsigned int thread=0; thread<n_threads; ++thread)
const Function<dim> &rhs,
Vector<double> &rhs_vector,
const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
+ const IteratorRange<DoFHandler<dim> > range,
Threads::ThreadMutex &mutex)
{
UpdateFlags update_flags = UpdateFlags(update_values |
}
-#if deal_II_dimension == 1
-void MatrixCreator::create_boundary_mass_matrix (const Mapping<1> &,
- const DoFHandler<1> &,
- const Quadrature<0> &,
- SparseMatrix<double> &,
- const FunctionMap<1>::type&,
- Vector<double> &,
- std::vector<unsigned int> &,
- const Function<1> * const)
+template <int dim, typename number>
+void MatrixCreator::create_mass_matrix (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<number> &matrix,
+ const Function<dim> * const coefficient)
{
- // what would that be in 1d? the
- // identity matrix on the boundary
- // dofs?
- Assert (false, ExcNotImplemented());
-}
-
-
-#endif
-
-
+ Assert (matrix.m() == dof.n_dofs(),
+ ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+ Assert (matrix.n() == dof.n_dofs(),
+ ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-template <int dim>
-void
-MatrixCreator::create_boundary_mass_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
- SparseMatrix<double> &matrix,
- const typename FunctionMap<dim>::type &boundary_functions,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const coefficient)
-{
const unsigned int n_threads = multithread_info.n_default_threads;
Threads::ThreadGroup<> threads;
// define starting and end point
// for each thread
- typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ typedef typename hp::DoFHandler<dim>::active_cell_iterator active_cell_iterator;
std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
= Threads::split_range<active_cell_iterator> (dof.begin_active(),
dof.end(), n_threads);
Threads::ThreadMutex mutex;
// then assemble in parallel
- typedef void (*create_boundary_mass_matrix_1_t)
- (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
- SparseMatrix<double> &matrix,
- const typename FunctionMap<dim>::type &boundary_functions,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
- Threads::ThreadMutex &mutex);
- create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1<dim>;
+ typedef void (*create_mass_matrix_1_t) (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<number> &matrix,
+ const Function<dim> * const coefficient,
+ const IteratorRange<hp::DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex);
+ create_mass_matrix_1_t p = &MatrixCreator::template create_mass_matrix_1<dim>;
for (unsigned int thread=0; thread<n_threads; ++thread)
- threads += Threads::spawn (p)(mapping, dof, q, matrix,
- boundary_functions, rhs_vector,
- dof_to_boundary_mapping, coefficient,
+ threads += Threads::spawn (p)(mapping, dof, q, matrix, coefficient,
thread_ranges[thread], mutex);
threads.join_all ();
}
-template <int dim>
+template <int dim, typename number>
void
-MatrixCreator::
-create_boundary_mass_matrix_1 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
- SparseMatrix<double> &matrix,
- const typename FunctionMap<dim>::type &boundary_functions,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
- Threads::ThreadMutex &mutex)
+MatrixCreator::create_mass_matrix_1 (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<number> &matrix,
+ const Function<dim> * const coefficient,
+ const IteratorRange<hp::DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex)
{
- const FiniteElement<dim> &fe = dof.get_fe();
- const unsigned int n_components = fe.n_components();
- const bool fe_is_system = (n_components != 1);
-
- 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.get_fe() == fe, ExcInternalError());
- Assert (dof_to_boundary_mapping.size() == dof.n_dofs(),
- ExcInternalError());
- Assert (n_components == boundary_functions.begin()->second->n_components,
- ExcComponentMismatch());
- Assert (coefficient ==0 ||
- coefficient->n_components==1 ||
- coefficient->n_components==n_components, ExcComponentMismatch());
-#ifdef DEBUG
- if (true)
- {
- unsigned int max_element = 0;
- for (std::vector<unsigned int>::const_iterator i=dof_to_boundary_mapping.begin();
- i!=dof_to_boundary_mapping.end(); ++i)
- if ((*i != DoFHandler<dim>::invalid_dof_index) &&
- (*i > max_element))
- max_element = *i;
- Assert (max_element == matrix.n()-1, ExcInternalError());
- };
-#endif
-
- const unsigned int dofs_per_cell = fe.dofs_per_cell,
- dofs_per_face = fe.dofs_per_face;
-
- FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
- Vector<double> cell_vector(dofs_per_cell);
-
+ UpdateFlags update_flags = UpdateFlags(update_values |
+ update_JxW_values);
+ if (coefficient != 0)
+ update_flags = UpdateFlags (update_flags | update_q_points);
- UpdateFlags update_flags = UpdateFlags (update_values |
- update_JxW_values |
- update_q_points);
- FEFaceValues<dim> fe_values (mapping, fe, q, update_flags);
+ hp::FEValues<dim> x_fe_values (mapping, dof.get_fe(), q, update_flags);
+
+ const unsigned int n_components = dof.get_fe().n_components();
- // two variables for the coefficient,
- // one for the two cases indicated in
- // the name
- std::vector<double> coefficient_values (fe_values.n_quadrature_points);
- std::vector<Vector<double> > coefficient_vector_values (fe_values.n_quadrature_points,
- Vector<double>(n_components));
+ Assert(coefficient == 0 ||
+ coefficient->n_components==1 ||
+ coefficient->n_components==n_components, ExcComponentMismatch());
- 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_components));
+ FullMatrix<double> cell_matrix;
+ std::vector<double> coefficient_values;
+ std::vector<Vector<double> > coefficient_vector_values;
+
+ std::vector<unsigned int> dof_indices;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator cell = range.first;
+ for (; cell!=range.second; ++cell)
+ {
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values ();
- std::vector<unsigned int> dofs (dofs_per_cell);
- std::vector<unsigned int> dofs_on_face_vector (dofs_per_face);
+ const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+ n_q_points = fe_values.n_quadrature_points;
+ const FiniteElement<dim> &fe = fe_values.get_fe();
- // for each dof on the cell, have a
- // flag whether it is on the face
- std::vector<bool> dof_is_on_face(dofs_per_cell);
+ cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+ coefficient_values.resize (n_q_points);
+ coefficient_vector_values.resize (n_q_points,
+ Vector<double> (n_components));
+ dof_indices.resize (dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator cell = range.first;
- for (; cell!=range.second; ++cell)
- 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_indicator()) !=
- boundary_functions.end())
+
+ cell_matrix = 0;
+ cell->get_dof_indices (dof_indices);
+
+ if (coefficient != 0)
{
- cell_matrix = 0;
- cell_vector = 0;
-
- fe_values.reinit (cell, face);
-
- if (fe_is_system)
- // FE has several components
+ if (coefficient->n_components==1)
{
- boundary_functions.find(cell->face(face)->boundary_indicator())
- ->second->vector_value_list (fe_values.get_quadrature_points(),
- rhs_values_system);
-
- if (coefficient != 0)
+ coefficient->value_list (fe_values.get_quadrature_points(),
+ coefficient_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
{
- if (coefficient->n_components==1)
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
{
- coefficient->value_list (fe_values.get_quadrature_points(),
- coefficient_values);
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
{
- 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);
- if (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first)
- {
- cell_matrix(i,j)
- += (u * v * weight * coefficient_values[point]);
- }
- }
- cell_vector(i) += v *
- rhs_values_system[point](
- fe.system_to_component_index(i).first) * weight;
- }
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (fe_values.shape_value(i,point) *
+ fe_values.shape_value(j,point) *
+ weight *
+ coefficient_values[point]);
}
}
- else
+ }
+ }
+ else
+ {
+ coefficient->vector_value_list (fe_values.get_quadrature_points(),
+ coefficient_vector_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
{
- 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 unsigned int component_i=
+ fe.system_to_component_index(i).first;
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
{
- 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)
- {
- const double u = fe_values.shape_value(j,point);
- if (fe.system_to_component_index(j).first ==
- component_i)
- {
- cell_matrix(i,j) +=
- (u * v * weight * coefficient_vector_values[point](component_i));
- }
- }
- cell_vector(i) += v * rhs_values_system[point](component_i) * weight;
- }
+ if ((n_components==1) ||
+ (fe.system_to_component_index(j).first == component_i))
+ cell_matrix(i,j) += (fe_values.shape_value(i,point) *
+ fe_values.shape_value(j,point) *
+ weight *
+ coefficient_vector_values[point](component_i));
}
}
}
- 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)
- {
- const double u = fe_values.shape_value(j,point);
- if (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first)
- {
- cell_matrix(i,j) += (u * v * weight);
- }
- }
- cell_vector(i) += v *
- rhs_values_system[point](
- fe.system_to_component_index(i).first) *
- weight;
- }
- }
}
- else
+ }
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (fe_values.shape_value(i,point) *
+ fe_values.shape_value(j,point) *
+ weight);
+ }
+ }
+ }
+
+ // transfer everything into the
+ // global object. lock the
+ // matrix meanwhile
+ Threads::ThreadMutex::ScopedLock lock (mutex);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ matrix.add (dof_indices[i], dof_indices[j],
+ cell_matrix(i,j));
+ }
+}
+
+
+
+template <int dim, typename number>
+void MatrixCreator::create_mass_matrix (const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<number> &matrix,
+ const Function<dim> * const coefficient)
+{
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ create_mass_matrix(hp::StaticMappingQ1<dim>::mapping_collection, dof, q, matrix, coefficient);
+}
+
+
+
+template <int dim, typename number>
+void MatrixCreator::create_mass_matrix (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<number> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient)
+{
+ Assert (matrix.m() == dof.n_dofs(),
+ ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+ Assert (matrix.n() == dof.n_dofs(),
+ ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ Threads::ThreadGroup<> threads;
+
+ // define starting and end point
+ // for each thread
+ typedef typename hp::DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
+ = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+ dof.end(), n_threads);
+
+ // mutex to synchronise access to
+ // the matrix
+ Threads::ThreadMutex mutex;
+
+ // then assemble in parallel
+ typedef void (*create_mass_matrix_2_t) (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<number> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient,
+ const IteratorRange<hp::DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex);
+ create_mass_matrix_2_t p = &MatrixCreator::template create_mass_matrix_2<dim>;
+ for (unsigned int thread=0; thread<n_threads; ++thread)
+ threads += Threads::spawn (p)(mapping, dof, q, matrix, rhs,
+ rhs_vector, coefficient,
+ thread_ranges[thread], mutex);
+ threads.join_all ();
+}
+
+
+
+template <int dim, typename number>
+void
+MatrixCreator::create_mass_matrix_2 (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<number> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient,
+ const IteratorRange<hp::DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex)
+{
+ UpdateFlags update_flags = UpdateFlags(update_values |
+ update_q_points |
+ update_JxW_values);
+ if (coefficient != 0)
+ update_flags = UpdateFlags (update_flags | update_q_points);
+
+ hp::FEValues<dim> x_fe_values (mapping, dof.get_fe(), q, update_flags);
+
+ const unsigned int n_components = dof.get_fe().n_components();
+
+ Assert(coefficient == 0 ||
+ coefficient->n_components==1 ||
+ coefficient->n_components==n_components, ExcComponentMismatch());
+
+ FullMatrix<double> cell_matrix;
+ Vector<double> local_rhs;
+ std::vector<double> rhs_values;
+ std::vector<double> coefficient_values;
+ std::vector<Vector<double> > coefficient_vector_values;
+
+ std::vector<unsigned int> dof_indices;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator cell = range.first;
+ for (; cell!=range.second; ++cell)
+ {
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = 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> &fe = fe_values.get_fe();
+
+ cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+ local_rhs.reinit (dofs_per_cell);
+ rhs_values.resize (fe_values.n_quadrature_points);
+ coefficient_values.resize (n_q_points);
+ coefficient_vector_values.resize (n_q_points,
+ Vector<double> (n_components));
+ dof_indices.resize (dofs_per_cell);
+
+
+ cell_matrix = 0;
+ local_rhs = 0;
+ cell->get_dof_indices (dof_indices);
+
+ rhs.value_list (fe_values.get_quadrature_points(), rhs_values);
+
+ if (coefficient != 0)
+ {
+ if (coefficient->n_components==1)
+ {
+ coefficient->value_list (fe_values.get_quadrature_points(),
+ coefficient_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (fe_values.shape_value(i,point) *
+ fe_values.shape_value(j,point) *
+ weight *
+ coefficient_values[point]);
+ local_rhs(i) += fe_values.shape_value(i,point) *
+ rhs_values[point] * weight;
+ }
+ }
+ }
+ }
+ else
+ {
+ coefficient->vector_value_list (fe_values.get_quadrature_points(),
+ coefficient_vector_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const unsigned int component_i=
+ fe.system_to_component_index(i).first;
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ if ((n_components==1) ||
+ (fe.system_to_component_index(j).first == component_i))
+ cell_matrix(i,j) += (fe_values.shape_value(i,point) *
+ fe_values.shape_value(j,point) *
+ weight *
+ coefficient_vector_values[point](component_i));
+ local_rhs(i) += fe_values.shape_value(i,point) *
+ rhs_values[point] * weight;
+ }
+ }
+ }
+ }
+ }
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const double v = fe_values.shape_value(i,point);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (fe_values.shape_value(i,point) *
+ fe_values.shape_value(j,point) *
+ weight);
+ local_rhs(i) += v * rhs_values[point] * weight;
+ }
+ }
+ }
+
+ // transfer everything into the
+ // global object. lock the
+ // matrix meanwhile
+ Threads::ThreadMutex::ScopedLock lock (mutex);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ matrix.add (dof_indices[i], dof_indices[j],
+ cell_matrix(i,j));
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ rhs_vector(dof_indices[i]) += local_rhs(i);
+ }
+}
+
+
+
+template <int dim, typename number>
+void MatrixCreator::create_mass_matrix (const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<number> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient)
+{
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ create_mass_matrix(hp::StaticMappingQ1<dim>::mapping_collection, dof, q,
+ matrix, rhs, rhs_vector, coefficient);
+}
+
+
+
+
+#if deal_II_dimension == 1
+
+void MatrixCreator::create_boundary_mass_matrix (const Mapping<1> &,
+ const DoFHandler<1> &,
+ const Quadrature<0> &,
+ SparseMatrix<double> &,
+ const FunctionMap<1>::type&,
+ Vector<double> &,
+ std::vector<unsigned int> &,
+ const Function<1> * const)
+{
+ // what would that be in 1d? the
+ // identity matrix on the boundary
+ // dofs?
+ Assert (false, ExcNotImplemented());
+}
+
+
+#endif
+
+
+
+template <int dim>
+void
+MatrixCreator::create_boundary_mass_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const coefficient)
+{
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ Threads::ThreadGroup<> threads;
+
+ // define starting and end point
+ // for each thread
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
+ = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+ dof.end(), n_threads);
+
+ // mutex to synchronise access to
+ // the matrix
+ Threads::ThreadMutex mutex;
+
+ // then assemble in parallel
+ typedef void (*create_boundary_mass_matrix_1_t)
+ (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const coefficient,
+ const IteratorRange<DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex);
+ create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1<dim>;
+ for (unsigned int thread=0; thread<n_threads; ++thread)
+ threads += Threads::spawn (p)(mapping, dof, q, matrix,
+ boundary_functions, rhs_vector,
+ dof_to_boundary_mapping, coefficient,
+ thread_ranges[thread], mutex);
+ threads.join_all ();
+}
+
+
+
+template <int dim>
+void
+MatrixCreator::
+create_boundary_mass_matrix_1 (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const coefficient,
+ const IteratorRange<DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex)
+{
+ const FiniteElement<dim> &fe = dof.get_fe();
+ const unsigned int n_components = fe.n_components();
+ const bool fe_is_system = (n_components != 1);
+
+ 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.get_fe() == fe, ExcInternalError());
+ Assert (dof_to_boundary_mapping.size() == dof.n_dofs(),
+ ExcInternalError());
+ Assert (n_components == boundary_functions.begin()->second->n_components,
+ ExcComponentMismatch());
+ Assert (coefficient ==0 ||
+ coefficient->n_components==1 ||
+ coefficient->n_components==n_components, ExcComponentMismatch());
+#ifdef DEBUG
+ if (true)
+ {
+ unsigned int max_element = 0;
+ for (std::vector<unsigned int>::const_iterator i=dof_to_boundary_mapping.begin();
+ i!=dof_to_boundary_mapping.end(); ++i)
+ if ((*i != DoFHandler<dim>::invalid_dof_index) &&
+ (*i > max_element))
+ max_element = *i;
+ Assert (max_element == matrix.n()-1, ExcInternalError());
+ };
+#endif
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell,
+ dofs_per_face = fe.dofs_per_face;
+
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_vector(dofs_per_cell);
+
+
+ UpdateFlags update_flags = UpdateFlags (update_values |
+ update_JxW_values |
+ update_q_points);
+ FEFaceValues<dim> 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);
+ std::vector<Vector<double> > coefficient_vector_values (fe_values.n_quadrature_points,
+ Vector<double>(n_components));
+
+ 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_components));
+
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ std::vector<unsigned int> dofs_on_face_vector (dofs_per_face);
+
+ // for each dof on the cell, have a
+ // flag whether it is on the face
+ std::vector<bool> dof_is_on_face(dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+ for (; cell!=range.second; ++cell)
+ 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_indicator()) !=
+ boundary_functions.end())
+ {
+ cell_matrix = 0;
+ cell_vector = 0;
+
+ fe_values.reinit (cell, face);
+
+ if (fe_is_system)
+ // FE has several components
+ {
+ boundary_functions.find(cell->face(face)->boundary_indicator())
+ ->second->vector_value_list (fe_values.get_quadrature_points(),
+ rhs_values_system);
+
+ if (coefficient != 0)
+ {
+ if (coefficient->n_components==1)
+ {
+ 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);
+ if (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first)
+ {
+ cell_matrix(i,j)
+ += (u * v * weight * coefficient_values[point]);
+ }
+ }
+ cell_vector(i) += v *
+ rhs_values_system[point](
+ fe.system_to_component_index(i).first) * weight;
+ }
+ }
+ }
+ else
+ {
+ 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)
+ {
+ const double u = fe_values.shape_value(j,point);
+ if (fe.system_to_component_index(j).first ==
+ component_i)
+ {
+ cell_matrix(i,j) +=
+ (u * v * weight * coefficient_vector_values[point](component_i));
+ }
+ }
+ cell_vector(i) += v * rhs_values_system[point](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)
+ {
+ const double u = fe_values.shape_value(j,point);
+ if (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first)
+ {
+ cell_matrix(i,j) += (u * v * weight);
+ }
+ }
+ cell_vector(i) += v *
+ rhs_values_system[point](
+ fe.system_to_component_index(i).first) *
+ weight;
+ }
+ }
+ }
+ else
+ // FE is a scalar one
+ {
+ boundary_functions.find(cell->face(face)->boundary_indicator())
+ ->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);
+ cell_matrix(i,j) += (u * v * weight * coefficient_values[point]);
+ }
+ cell_vector(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);
+ cell_matrix(i,j) += (u * v * weight);
+ }
+ cell_vector(i) += v * rhs_values_scalar[point] * weight;
+ }
+ }
+ }
+
+ // 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.
+ cell->get_dof_indices (dofs);
+ cell->face(face)->get_dof_indices (dofs_on_face_vector);
+
+ // check for each of the
+ // dofs on this cell
+ // whether it is on the
+ // face
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ dof_is_on_face[i] = (std::find(dofs_on_face_vector.begin(),
+ dofs_on_face_vector.end(),
+ dofs[i])
+ !=
+ dofs_on_face_vector.end());
+
+ // 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.
+#ifdef DEBUG
+ double max_diag_entry = 0;
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (std::fabs(cell_matrix(i,i)) > max_diag_entry)
+ max_diag_entry = std::fabs(cell_matrix(i,i));
+#endif
+
+ // lock the matrix
+ Threads::ThreadMutex::ScopedLock lock (mutex);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (dof_is_on_face[i] && dof_is_on_face[j])
+ matrix.add(dof_to_boundary_mapping[dofs[i]],
+ dof_to_boundary_mapping[dofs[j]],
+ cell_matrix(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(cell_matrix(i,j)) <= 1e-10 * max_diag_entry,
+ ExcInternalError ());
+ };
+
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (dof_is_on_face[j])
+ rhs_vector(dof_to_boundary_mapping[dofs[j]]) += cell_vector(j);
+ else
+ {
+ // compare here for relative
+ // smallness
+ Assert (std::fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
+ ExcInternalError());
+ };
+ };
+}
+
+
+template <int dim>
+void MatrixCreator::create_boundary_mass_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &rhs,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const a)
+{
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ create_boundary_mass_matrix(StaticMappingQ1<dim>::mapping, dof, q,
+ matrix,rhs, rhs_vector, dof_to_boundary_mapping, a);
+}
+
+
+
+#if deal_II_dimension == 1
+
+void MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection<1> &,
+ const hp::DoFHandler<1> &,
+ const hp::QCollection<0> &,
+ SparseMatrix<double> &,
+ const FunctionMap<1>::type&,
+ Vector<double> &,
+ std::vector<unsigned int> &,
+ const Function<1> * const)
+{
+ // what would that be in 1d? the
+ // identity matrix on the boundary
+ // dofs?
+ Assert (false, ExcNotImplemented());
+}
+
+
+#endif
+
+
+
+template <int dim>
+void
+MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const coefficient)
+{
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ Threads::ThreadGroup<> threads;
+
+ // define starting and end point
+ // for each thread
+ typedef typename hp::DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
+ = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+ dof.end(), n_threads);
+
+ // mutex to synchronise access to
+ // the matrix
+ Threads::ThreadMutex mutex;
+
+ // then assemble in parallel
+ typedef void (*create_boundary_mass_matrix_1_t)
+ (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const coefficient,
+ const IteratorRange<hp::DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex);
+ create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1<dim>;
+ for (unsigned int thread=0; thread<n_threads; ++thread)
+ threads += Threads::spawn (p)(mapping, dof, q, matrix,
+ boundary_functions, rhs_vector,
+ dof_to_boundary_mapping, coefficient,
+ thread_ranges[thread], mutex);
+ threads.join_all ();
+}
+
+
+
+template <int dim>
+void
+MatrixCreator::
+create_boundary_mass_matrix_1 (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const coefficient,
+ const IteratorRange<hp::DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex)
+{
+/*
+ const FiniteElement<dim> &fe = dof.get_fe();
+ const unsigned int n_components = fe.n_components();
+ const bool fe_is_system = (n_components != 1);
+
+ 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.get_fe() == fe, ExcInternalError());
+ Assert (dof_to_boundary_mapping.size() == dof.n_dofs(),
+ ExcInternalError());
+ Assert (n_components == boundary_functions.begin()->second->n_components,
+ ExcComponentMismatch());
+ Assert (coefficient ==0 ||
+ coefficient->n_components==1 ||
+ coefficient->n_components==n_components, ExcComponentMismatch());
+#ifdef DEBUG
+ if (true)
+ {
+ unsigned int max_element = 0;
+ for (std::vector<unsigned int>::const_iterator i=dof_to_boundary_mapping.begin();
+ i!=dof_to_boundary_mapping.end(); ++i)
+ if ((*i != hp::DoFHandler<dim>::invalid_dof_index) &&
+ (*i > max_element))
+ max_element = *i;
+ Assert (max_element == matrix.n()-1, ExcInternalError());
+ };
+#endif
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell,
+ dofs_per_face = fe.dofs_per_face;
+
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_vector(dofs_per_cell);
+
+
+ UpdateFlags update_flags = UpdateFlags (update_values |
+ update_JxW_values |
+ update_q_points);
+ FEFaceValues<dim> 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);
+ std::vector<Vector<double> > coefficient_vector_values (fe_values.n_quadrature_points,
+ Vector<double>(n_components));
+
+ 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_components));
+
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ std::vector<unsigned int> dofs_on_face_vector (dofs_per_face);
+
+ // for each dof on the cell, have a
+ // flag whether it is on the face
+ std::vector<bool> dof_is_on_face(dofs_per_cell);
+
+ typename hp::DoFHandler<dim>::active_cell_iterator cell = range.first;
+ for (; cell!=range.second; ++cell)
+ 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_indicator()) !=
+ boundary_functions.end())
+ {
+ cell_matrix = 0;
+ cell_vector = 0;
+
+ fe_values.reinit (cell, face);
+
+ if (fe_is_system)
+ // FE has several components
+ {
+ boundary_functions.find(cell->face(face)->boundary_indicator())
+ ->second->vector_value_list (fe_values.get_quadrature_points(),
+ rhs_values_system);
+
+ if (coefficient != 0)
+ {
+ if (coefficient->n_components==1)
+ {
+ 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);
+ if (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first)
+ {
+ cell_matrix(i,j)
+ += (u * v * weight * coefficient_values[point]);
+ }
+ }
+ cell_vector(i) += v *
+ rhs_values_system[point](
+ fe.system_to_component_index(i).first) * weight;
+ }
+ }
+ }
+ else
+ {
+ 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)
+ {
+ const double u = fe_values.shape_value(j,point);
+ if (fe.system_to_component_index(j).first ==
+ component_i)
+ {
+ cell_matrix(i,j) +=
+ (u * v * weight * coefficient_vector_values[point](component_i));
+ }
+ }
+ cell_vector(i) += v * rhs_values_system[point](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)
+ {
+ const double u = fe_values.shape_value(j,point);
+ if (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first)
+ {
+ cell_matrix(i,j) += (u * v * weight);
+ }
+ }
+ cell_vector(i) += v *
+ rhs_values_system[point](
+ fe.system_to_component_index(i).first) *
+ weight;
+ }
+ }
+ }
+ else
// FE is a scalar one
{
- boundary_functions.find(cell->face(face)->boundary_indicator())
- ->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);
- cell_matrix(i,j) += (u * v * weight * coefficient_values[point]);
- }
- cell_vector(i) += v * rhs_values_scalar[point] *weight;
- }
- }
- }
- else
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ boundary_functions.find(cell->face(face)->boundary_indicator())
+ ->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);
+ cell_matrix(i,j) += (u * v * weight * coefficient_values[point]);
+ }
+ cell_vector(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);
+ cell_matrix(i,j) += (u * v * weight);
+ }
+ cell_vector(i) += v * rhs_values_scalar[point] * weight;
+ }
+ }
+ }
+
+ // 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.
+ cell->get_dof_indices (dofs);
+ cell->face(face)->get_dof_indices (dofs_on_face_vector);
+
+ // check for each of the
+ // dofs on this cell
+ // whether it is on the
+ // face
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ dof_is_on_face[i] = (std::find(dofs_on_face_vector.begin(),
+ dofs_on_face_vector.end(),
+ dofs[i])
+ !=
+ dofs_on_face_vector.end());
+
+ // 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.
+#ifdef DEBUG
+ double max_diag_entry = 0;
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (std::fabs(cell_matrix(i,i)) > max_diag_entry)
+ max_diag_entry = std::fabs(cell_matrix(i,i));
+#endif
+
+ // lock the matrix
+ Threads::ThreadMutex::ScopedLock lock (mutex);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (dof_is_on_face[i] && dof_is_on_face[j])
+ matrix.add(dof_to_boundary_mapping[dofs[i]],
+ dof_to_boundary_mapping[dofs[j]],
+ cell_matrix(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(cell_matrix(i,j)) <= 1e-10 * max_diag_entry,
+ ExcInternalError ());
+ };
+
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (dof_is_on_face[j])
+ rhs_vector(dof_to_boundary_mapping[dofs[j]]) += cell_vector(j);
+ else
+ {
+ // compare here for relative
+ // smallness
+ Assert (std::fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
+ ExcInternalError());
+ };
+ };
+*/
+}
+
+
+template <int dim>
+void MatrixCreator::create_boundary_mass_matrix (const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &rhs,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const a)
+{
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ create_boundary_mass_matrix(hp::StaticMappingQ1<dim>::mapping_collection, dof, q,
+ matrix,rhs, rhs_vector, dof_to_boundary_mapping, a);
+}
+
+
+
+template <int dim>
+void MatrixCreator::create_laplace_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient)
+{
+ Assert (matrix.m() == dof.n_dofs(),
+ ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+ Assert (matrix.n() == dof.n_dofs(),
+ ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ Threads::ThreadGroup<> threads;
+
+ // define starting and end point
+ // for each thread
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
+ = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+ dof.end(), n_threads);
+
+ // mutex to synchronise access to
+ // the matrix
+ Threads::ThreadMutex mutex;
+
+ // then assemble in parallel
+ typedef void (*create_laplace_matrix_1_t) (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient,
+ const IteratorRange<DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex);
+ create_laplace_matrix_1_t p = &MatrixCreator::template create_laplace_matrix_1<dim>;
+ for (unsigned int thread=0; thread<n_threads; ++thread)
+ threads += Threads::spawn (p)(mapping, dof, q, matrix, coefficient,
+ thread_ranges[thread], mutex);
+ threads.join_all ();
+}
+
+
+
+template <int dim>
+void MatrixCreator::create_laplace_matrix_1 (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient,
+ const IteratorRange<DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex)
+{
+ UpdateFlags update_flags = UpdateFlags(update_JxW_values |
+ update_gradients);
+ if (coefficient != 0)
+ update_flags = UpdateFlags (update_flags | update_q_points);
+
+ FEValues<dim> fe_values (mapping, dof.get_fe(), q, update_flags);
+
+ const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+ n_q_points = fe_values.n_quadrature_points;
+ const FiniteElement<dim> &fe = fe_values.get_fe();
+ const unsigned int n_components = fe.n_components();
+
+ Assert(coefficient == 0 ||
+ coefficient->n_components==1 ||
+ coefficient->n_components==n_components, ExcComponentMismatch());
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ std::vector<double> coefficient_values (n_q_points);
+ std::vector<Vector<double> > coefficient_vector_values (n_q_points,
+ Vector<double> (n_components));
+
+ std::vector<unsigned int> dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+ for (; cell!=range.second; ++cell)
+ {
+ fe_values.reinit (cell);
+
+ cell_matrix = 0;
+ cell->get_dof_indices (dof_indices);
+
+ if (coefficient != 0)
+ {
+ if (coefficient->n_components==1)
+ {
+ coefficient->value_list (fe_values.get_quadrature_points(),
+ coefficient_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (Du * Dv * weight *
+ coefficient_values[point]);
+ }
+ }
+ }
+ }
+ else
+ {
+ coefficient->vector_value_list (fe_values.get_quadrature_points(),
+ coefficient_vector_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
+ const unsigned int component_i=
+ fe.system_to_component_index(i).first;
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(j).first == component_i))
+ cell_matrix(i,j) += (Du * Dv * weight *
+ coefficient_vector_values[point](component_i));
+
+ }
+ }
+ }
+ }
+ }
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (Du * Dv * weight);
+ }
+ }
+ }
+
+ // transfer everything into the
+ // global object. lock the
+ // matrix meanwhile
+ Threads::ThreadMutex::ScopedLock lock (mutex);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ matrix.add (dof_indices[i], dof_indices[j],
+ cell_matrix(i,j));
+ }
+}
+
+
+
+template <int dim>
+void MatrixCreator::create_laplace_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient)
+{
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ create_laplace_matrix(StaticMappingQ1<dim>::mapping, dof, q, matrix, coefficient);
+}
+
+
+
+template <int dim>
+void MatrixCreator::create_laplace_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient)
+{
+ Assert (matrix.m() == dof.n_dofs(),
+ ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+ Assert (matrix.n() == dof.n_dofs(),
+ ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ Threads::ThreadGroup<> threads;
+
+ // define starting and end point
+ // for each thread
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
+ = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+ dof.end(), n_threads);
+
+ // mutex to synchronise access to
+ // the matrix
+ Threads::ThreadMutex mutex;
+
+ // then assemble in parallel
+ typedef void (*create_laplace_matrix_2_t) (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient,
+ const IteratorRange<DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex);
+ create_laplace_matrix_2_t p = &MatrixCreator::template create_laplace_matrix_2<dim>;
+ for (unsigned int thread=0; thread<n_threads; ++thread)
+ threads += Threads::spawn (p)(mapping, dof, q, matrix, rhs,
+ rhs_vector, coefficient,
+ thread_ranges[thread], mutex);
+ threads.join_all ();
+}
+
+
+
+template <int dim>
+void
+MatrixCreator::create_laplace_matrix_2 (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient,
+ const IteratorRange<DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex)
+{
+ UpdateFlags update_flags = UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values);
+ if (coefficient != 0)
+ update_flags = UpdateFlags (update_flags | update_q_points);
+
+ FEValues<dim> fe_values (mapping, dof.get_fe(), q, update_flags);
+
+ const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+ n_q_points = fe_values.n_quadrature_points;
+ const FiniteElement<dim> &fe = fe_values.get_fe();
+ const unsigned int n_components = fe.n_components();
+
+ Assert(coefficient == 0 ||
+ coefficient->n_components==1 ||
+ coefficient->n_components==n_components, ExcComponentMismatch());
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> local_rhs (dofs_per_cell);
+ std::vector<double> rhs_values (fe_values.n_quadrature_points);
+ std::vector<double> coefficient_values (n_q_points);
+ std::vector<Vector<double> > coefficient_vector_values (n_q_points,
+ Vector<double> (n_components));
+
+ std::vector<unsigned int> dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+ for (; cell!=range.second; ++cell)
+ {
+ fe_values.reinit (cell);
+
+ cell_matrix = 0;
+ local_rhs = 0;
+ cell->get_dof_indices (dof_indices);
+
+ rhs.value_list (fe_values.get_quadrature_points(), rhs_values);
+
+ if (coefficient != 0)
+ {
+ if (coefficient->n_components==1)
+ {
+ coefficient->value_list (fe_values.get_quadrature_points(),
+ coefficient_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
{
const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+ for (unsigned int i=0; i<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 Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
{
- const double u = fe_values.shape_value(j,point);
- cell_matrix(i,j) += (u * v * weight);
+ const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (Du * Dv * weight *
+ coefficient_values[point]);
+ local_rhs(i) += v * rhs_values[point] * weight;
}
- cell_vector(i) += v * rhs_values_scalar[point] * weight;
}
}
}
-
- // 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.
- cell->get_dof_indices (dofs);
- cell->face(face)->get_dof_indices (dofs_on_face_vector);
-
- // check for each of the
- // dofs on this cell
- // whether it is on the
- // face
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- dof_is_on_face[i] = (std::find(dofs_on_face_vector.begin(),
- dofs_on_face_vector.end(),
- dofs[i])
- !=
- dofs_on_face_vector.end());
-
- // 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.
-#ifdef DEBUG
- double max_diag_entry = 0;
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- if (std::fabs(cell_matrix(i,i)) > max_diag_entry)
- max_diag_entry = std::fabs(cell_matrix(i,i));
-#endif
-
- // lock the matrix
- Threads::ThreadMutex::ScopedLock lock (mutex);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (dof_is_on_face[i] && dof_is_on_face[j])
- matrix.add(dof_to_boundary_mapping[dofs[i]],
- dof_to_boundary_mapping[dofs[j]],
- cell_matrix(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(cell_matrix(i,j)) <= 1e-10 * max_diag_entry,
- ExcInternalError ());
- };
-
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (dof_is_on_face[j])
- rhs_vector(dof_to_boundary_mapping[dofs[j]]) += cell_vector(j);
- else
+ else
+ {
+ coefficient->vector_value_list (fe_values.get_quadrature_points(),
+ coefficient_vector_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const double v = fe_values.shape_value(i,point);
+ const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
+ const unsigned int component_i=
+ fe.system_to_component_index(i).first;
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(j).first == component_i))
+ cell_matrix(i,j) += (Du * Dv * weight *
+ coefficient_vector_values[point](component_i));
+ local_rhs(i) += v * rhs_values[point] * weight;
+ }
+ }
+ }
+ }
+ }
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
{
- // compare here for relative
- // smallness
- Assert (std::fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
- ExcInternalError());
- };
- };
+ const double v = fe_values.shape_value(i,point);
+ const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (Du * Dv * weight);
+ local_rhs(i) += v * rhs_values[point] * weight;
+ }
+ }
+ }
+
+ // transfer everything into the
+ // global object. lock the
+ // matrix meanwhile
+ Threads::ThreadMutex::ScopedLock lock (mutex);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ matrix.add (dof_indices[i], dof_indices[j],
+ cell_matrix(i,j));
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ rhs_vector(dof_indices[i]) += local_rhs(i);
+ };
}
+
template <int dim>
-void MatrixCreator::create_boundary_mass_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
- SparseMatrix<double> &matrix,
- const typename FunctionMap<dim>::type &rhs,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const a)
+void MatrixCreator::create_laplace_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- create_boundary_mass_matrix(StaticMappingQ1<dim>::mapping, dof, q,
- matrix,rhs, rhs_vector, dof_to_boundary_mapping, a);
+ create_laplace_matrix(StaticMappingQ1<dim>::mapping, dof, q,
+ matrix, rhs, rhs_vector, coefficient);
}
template <int dim>
-void MatrixCreator::create_laplace_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
+void MatrixCreator::create_laplace_matrix (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> * const coefficient)
{
// define starting and end point
// for each thread
- typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ typedef typename hp::DoFHandler<dim>::active_cell_iterator active_cell_iterator;
std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
= Threads::split_range<active_cell_iterator> (dof.begin_active(),
dof.end(), n_threads);
Threads::ThreadMutex mutex;
// then assemble in parallel
- typedef void (*create_laplace_matrix_1_t) (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
+ typedef void (*create_laplace_matrix_1_t) (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
+ const IteratorRange<hp::DoFHandler<dim> > range,
Threads::ThreadMutex &mutex);
create_laplace_matrix_1_t p = &MatrixCreator::template create_laplace_matrix_1<dim>;
for (unsigned int thread=0; thread<n_threads; ++thread)
template <int dim>
-void MatrixCreator::create_laplace_matrix_1 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
- Threads::ThreadMutex &mutex)
+void
+MatrixCreator::create_laplace_matrix_1 (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient,
+ const IteratorRange<hp::DoFHandler<dim> > range,
+ Threads::ThreadMutex &mutex)
{
- UpdateFlags update_flags = UpdateFlags(update_JxW_values |
- update_gradients);
+ UpdateFlags update_flags = UpdateFlags(update_gradients |
+ update_JxW_values);
if (coefficient != 0)
update_flags = UpdateFlags (update_flags | update_q_points);
- FEValues<dim> fe_values (mapping, dof.get_fe(), q, update_flags);
+ hp::FEValues<dim> x_fe_values (mapping, dof.get_fe(), q, update_flags);
- const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
- n_q_points = fe_values.n_quadrature_points;
- const FiniteElement<dim> &fe = fe_values.get_fe();
- const unsigned int n_components = fe.n_components();
+ const unsigned int n_components = dof.get_fe().n_components();
Assert(coefficient == 0 ||
coefficient->n_components==1 ||
coefficient->n_components==n_components, ExcComponentMismatch());
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- std::vector<double> coefficient_values (n_q_points);
- std::vector<Vector<double> > coefficient_vector_values (n_q_points,
- Vector<double> (n_components));
+ FullMatrix<double> cell_matrix;
+ std::vector<double> coefficient_values;
+ std::vector<Vector<double> > coefficient_vector_values;
- std::vector<unsigned int> dof_indices (dofs_per_cell);
+ std::vector<unsigned int> dof_indices;
- typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+ typename hp::DoFHandler<dim>::active_cell_iterator cell = range.first;
for (; cell!=range.second; ++cell)
{
- fe_values.reinit (cell);
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = 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> &fe = fe_values.get_fe();
+
+ cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+ coefficient_values.resize (n_q_points);
+ coefficient_vector_values.resize (n_q_points,
+ Vector<double> (n_components));
+ dof_indices.resize (dofs_per_cell);
+
cell_matrix = 0;
cell->get_dof_indices (dof_indices);
coefficient->value_list (fe_values.get_quadrature_points(),
coefficient_values);
for (unsigned int point=0; point<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
- if ((n_components==1) ||
- (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first))
- cell_matrix(i,j) += (Du * Dv * weight *
- coefficient_values[point]);
- }
- }
- }
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) += (Du * Dv * weight *
+ coefficient_values[point]);
+ }
+ }
+ }
}
else
{
coefficient->vector_value_list (fe_values.get_quadrature_points(),
coefficient_vector_values);
for (unsigned int point=0; point<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
- const unsigned int component_i=
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first == component_i))
- cell_matrix(i,j) += (Du * Dv * weight *
- coefficient_vector_values[point](component_i));
-
- }
- }
- }
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
+ const unsigned int component_i=
+ fe.system_to_component_index(i).first;
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const Tensor<1,dim>& Du = fe_values.shape_grad(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(j).first == component_i))
+ cell_matrix(i,j) += (Du * Dv * weight *
+ coefficient_vector_values[point](component_i));
+ }
+ }
+ }
}
}
else
for (unsigned int point=0; point<n_q_points; ++point)
{
const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim>& Dv = fe_values.shape_grad(i,point);
for (unsigned int j=0; j<dofs_per_cell; ++j)
}
}
}
-
+
// transfer everything into the
// global object. lock the
// matrix meanwhile
template <int dim>
-void MatrixCreator::create_laplace_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
+void MatrixCreator::create_laplace_matrix (const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> * const coefficient)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- create_laplace_matrix(StaticMappingQ1<dim>::mapping, dof, q, matrix, coefficient);
+ create_laplace_matrix(hp::StaticMappingQ1<dim>::mapping_collection, dof, q, matrix, coefficient);
}
template <int dim>
-void MatrixCreator::create_laplace_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
+void MatrixCreator::create_laplace_matrix (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> &rhs,
Vector<double> &rhs_vector,
// define starting and end point
// for each thread
- typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ typedef typename hp::DoFHandler<dim>::active_cell_iterator active_cell_iterator;
std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
= Threads::split_range<active_cell_iterator> (dof.begin_active(),
dof.end(), n_threads);
Threads::ThreadMutex mutex;
// then assemble in parallel
- typedef void (*create_laplace_matrix_2_t) (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
+ typedef void (*create_laplace_matrix_2_t) (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> &rhs,
Vector<double> &rhs_vector,
const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
+ const IteratorRange<hp::DoFHandler<dim> > range,
Threads::ThreadMutex &mutex);
create_laplace_matrix_2_t p = &MatrixCreator::template create_laplace_matrix_2<dim>;
for (unsigned int thread=0; thread<n_threads; ++thread)
template <int dim>
void
-MatrixCreator::create_laplace_matrix_2 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
+MatrixCreator::create_laplace_matrix_2 (const hp::MappingCollection<dim> &mapping,
+ const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> &rhs,
Vector<double> &rhs_vector,
const Function<dim> * const coefficient,
- const IteratorRange<dim> range,
+ const IteratorRange<hp::DoFHandler<dim> > range,
Threads::ThreadMutex &mutex)
{
UpdateFlags update_flags = UpdateFlags(update_values |
if (coefficient != 0)
update_flags = UpdateFlags (update_flags | update_q_points);
- FEValues<dim> fe_values (mapping, dof.get_fe(), q, update_flags);
+ hp::FEValues<dim> x_fe_values (mapping, dof.get_fe(), q, update_flags);
- const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
- n_q_points = fe_values.n_quadrature_points;
- const FiniteElement<dim> &fe = fe_values.get_fe();
- const unsigned int n_components = fe.n_components();
+ const unsigned int n_components = dof.get_fe().n_components();
Assert(coefficient == 0 ||
coefficient->n_components==1 ||
coefficient->n_components==n_components, ExcComponentMismatch());
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> local_rhs (dofs_per_cell);
- std::vector<double> rhs_values (fe_values.n_quadrature_points);
- std::vector<double> coefficient_values (n_q_points);
- std::vector<Vector<double> > coefficient_vector_values (n_q_points,
- Vector<double> (n_components));
+ FullMatrix<double> cell_matrix;
+ Vector<double> local_rhs;
+ std::vector<double> rhs_values;
+ std::vector<double> coefficient_values;
+ std::vector<Vector<double> > coefficient_vector_values;
- std::vector<unsigned int> dof_indices (dofs_per_cell);
+ std::vector<unsigned int> dof_indices;
- typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+ typename hp::DoFHandler<dim>::active_cell_iterator cell = range.first;
for (; cell!=range.second; ++cell)
{
- fe_values.reinit (cell);
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = 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> &fe = fe_values.get_fe();
+
+ cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+ local_rhs.reinit (dofs_per_cell);
+ rhs_values.resize (fe_values.n_quadrature_points);
+ coefficient_values.resize (n_q_points);
+ coefficient_vector_values.resize (n_q_points,
+ Vector<double> (n_components));
+ dof_indices.resize (dofs_per_cell);
+
cell_matrix = 0;
local_rhs = 0;
cell_matrix(i,j));
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dof_indices[i]) += local_rhs(i);
- };
+ }
}
template <int dim>
-void MatrixCreator::create_laplace_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
+void MatrixCreator::create_laplace_matrix (const hp::DoFHandler<dim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> &rhs,
Vector<double> &rhs_vector,
const Function<dim> * const coefficient)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- create_laplace_matrix(StaticMappingQ1<dim>::mapping, dof, q,
+ create_laplace_matrix(hp::StaticMappingQ1<dim>::mapping_collection, dof, q,
matrix, rhs, rhs_vector, coefficient);
}
// explicit instantiations
+// non-hp version of create_mass_matrix
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
(const Mapping<deal_II_dimension> &mapping,
Vector<double> &rhs_vector,
const Function<deal_II_dimension> * const coefficient);
+// hp versions of functions
+template
+void MatrixCreator::create_mass_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);
+template
+void MatrixCreator::create_mass_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);
+template
+void MatrixCreator::create_mass_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);
+template
+void MatrixCreator::create_mass_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);
+
+
+template
+void MatrixCreator::create_mass_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<float> &matrix,
+ const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_matrix<deal_II_dimension>
+(const hp::DoFHandler<deal_II_dimension> &dof,
+ const hp::QCollection<deal_II_dimension> &q,
+ SparseMatrix<float> &matrix,
+ const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_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<float> &matrix,
+ const Function<deal_II_dimension> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_matrix<deal_II_dimension>
+(const hp::DoFHandler<deal_II_dimension> &dof,
+ const hp::QCollection<deal_II_dimension> &q,
+ SparseMatrix<float> &matrix,
+ const Function<deal_II_dimension> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<deal_II_dimension> * const coefficient);
+
+
#if deal_II_dimension != 1
template
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
const Function<deal_II_dimension> * const a);
+
+
+#if deal_II_dimension != 1
+template
+void
+MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
+(const hp::MappingCollection<deal_II_dimension> &mapping,
+ const hp::DoFHandler<deal_II_dimension> &dof,
+ const hp::QCollection<deal_II_dimension-1> &q,
+ SparseMatrix<double> &matrix,
+ const FunctionMap<deal_II_dimension>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<deal_II_dimension> * const a);
+#endif
+
+template
+void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
+(const hp::DoFHandler<deal_II_dimension> &dof,
+ const hp::QCollection<deal_II_dimension-1> &q,
+ SparseMatrix<double> &matrix,
+ const FunctionMap<deal_II_dimension>::type &rhs,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<deal_II_dimension> * const a);
+
+
+
+// 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);
-
template
void MatrixCreator::create_laplace_matrix<deal_II_dimension>
(const Mapping<deal_II_dimension> &mapping,
Vector<double> &rhs_vector,
const Function<deal_II_dimension> * const coefficient);
+// 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);
+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);
+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);
+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);
+
DEAL_II_NAMESPACE_CLOSE