#include <deal.II/base/function.h>
#include <deal.II/base/quadrature.h>
-#include <deal.II/base/thread_management.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/geometry_info.h>
data.dof_indices,
*matrix);
}
+
+
+
+ namespace AssemblerBoundary
+ {
+ struct Scratch
+ {
+ Scratch() {}
+ };
+
+ template <int dim,int spacedim>
+ 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 DoFHandler<dim,spacedim>::active_cell_iterator cell;
+ std::vector<FullMatrix<double> > cell_matrix;
+ std::vector<Vector<double> > cell_vector;
+ };
+
+ template <int dim,int spacedim>
+ CopyData<dim,spacedim>::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)
+ {}
+ }
}
}
{
template <int dim, int spacedim>
void
- create_boundary_mass_matrix_1 (std_cxx1x::tuple<const Mapping<dim, spacedim> &,
- const DoFHandler<dim,spacedim> &,
- const Quadrature<dim-1> & > commons,
- 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,
- const std::vector<unsigned int> &component_mapping,
- const MatrixCreator::internal::IteratorRange<DoFHandler<dim,spacedim> > range,
- Threads::Mutex &mutex)
+ create_boundary_mass_matrix_1 (typename DoFHandler<dim,spacedim>::active_cell_iterator const &cell,
+ MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<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 Mapping<dim,spacedim> &mapping = std_cxx1x::get<0>(commons);
- const DoFHandler<dim,spacedim> &dof = std_cxx1x::get<1>(commons);
- const Quadrature<dim-1>& q = std_cxx1x::get<2>(commons);
-
- const FiniteElement<dim,spacedim> &fe = dof.get_fe();
- const unsigned int n_components = fe.n_components();
+ 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_system = (n_components != 1);
const bool fe_is_primitive = fe.is_primitive();
- 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);
+ 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 |
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));
+ 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 (dofs_per_cell);
std::vector<types::global_dof_index> 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);
+ // 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_indicator()) !=
+ 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);
- typename DoFHandler<dim,spacedim>::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;
+ 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_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];
+ }
- fe_values.reinit (cell, face);
+ // 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.));
- 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);
+ 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 (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
- // geiven, 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
+ if (!base.conforms(FiniteElementData<dim>::H1) &&
+ base.conforms(FiniteElementData<dim>::Hdiv))
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)
- {
- cell_matrix(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);
- }
- }
- cell_vector(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)
- cell_matrix(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);
- cell_vector(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_indicator())
- ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
+ normal_adjustment[point][comp] = fe_values.normal_vector(point)(bcomp)
+ * fe_values.normal_vector(point)(bcomp);
+ }
- 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)
+ 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)
{
- 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]);
+ 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);
+ }
}
- cell_vector(i) += v * rhs_values_scalar[point] *weight;
+ 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;
}
- }
- }
- // 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);
- 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());
+ 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_indicator())
+ ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
- // lock the matrix
- Threads::Mutex::ScopedLock lock (mutex);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- if (dof_is_on_face[i] && dof_to_boundary_mapping[dofs[i]] != numbers::invalid_dof_index)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (dof_is_on_face[j] && dof_to_boundary_mapping[dofs[j]] != numbers::invalid_dof_index)
+ 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)
{
- Assert(numbers::is_finite(cell_matrix(i,j)), ExcNumberNotFinite());
- matrix.add(dof_to_boundary_mapping[dofs[i]],
- dof_to_boundary_mapping[dofs[j]],
- cell_matrix(i,j));
+ const double u = fe_values.shape_value(j,point);
+ copy_data.cell_matrix.back()(i,j) += (u*v*weight*coefficient_values[point]);
}
- Assert(numbers::is_finite(cell_vector(i)), ExcNumberNotFinite());
- rhs_vector(dof_to_boundary_mapping[dofs[i]]) += cell_vector(i);
- }
- }
- }
+ 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));
+ 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<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_indicator()) !=
+ 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)
+ {
+ Assert(numbers::is_finite(copy_data.cell_matrix[pos](i,j)),
+ ExcNumberNotFinite());
+ 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));
+ }
+ Assert(numbers::is_finite(copy_data.cell_vector[pos](i)), ExcNumberNotFinite());
+ rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += copy_data.cell_vector[pos](i);
+ }
+ }
+ ++pos;
+ }
+ }
}
+
template <>
void
- create_boundary_mass_matrix_1<2,3> (std_cxx1x::tuple<const Mapping<2,3> &,
- const DoFHandler<2,3> &,
- const Quadrature<1> & > ,
- SparseMatrix<double> &,
- const FunctionMap<3>::type &,
- Vector<double> &,
- std::vector<types::global_dof_index> &,
- const Function<3> *const ,
- const std::vector<unsigned int> &,
- const MatrixCreator::internal::IteratorRange<DoFHandler<2,3> > ,
- Threads::Mutex &)
+ create_boundary_mass_matrix_1<2,3> (typename DoFHandler<2,3>::active_cell_iterator const &cell,
+ MatrixCreator::internal::AssemblerBoundary::Scratch const
+ &scratch,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<2,3>
+ ©_data,
+ Mapping<2,3> const &mapping,
+ FiniteElement<2,3> const &fe,
+ Quadrature<1> const &q,
+ typename FunctionMap<3>::type const &boundary_functions,
+ Function<3> const *const coefficient,
+ std::vector<unsigned int> const &component_mapping)
{
Assert(false,ExcNotImplemented());
}
template <>
void
- create_boundary_mass_matrix_1<1,3> (std_cxx1x::tuple<const Mapping<1,3> &,
- const DoFHandler<1,3> &,
- const Quadrature<0> & > ,
- SparseMatrix<double> &,
- const FunctionMap<3>::type &,
- Vector<double> &,
- std::vector<types::global_dof_index> &,
- const Function<3> *const ,
- const std::vector<unsigned int> &,
- const MatrixCreator::internal::IteratorRange<DoFHandler<1,3> > ,
- Threads::Mutex &)
+ create_boundary_mass_matrix_1<1,3> (typename DoFHandler<1,3>::active_cell_iterator const &cell,
+ MatrixCreator::internal::AssemblerBoundary::Scratch const
+ &scratch,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<1,3>
+ ©_data,
+ Mapping<1,3> const &mapping,
+ FiniteElement<1,3> const &fe,
+ Quadrature<0> const &q,
+ typename FunctionMap<3>::type const &boundary_functions,
+ Function<3> const *const coefficient,
+ std::vector<unsigned int> const &component_mapping)
{
Assert(false,ExcNotImplemented());
}
else
AssertDimension (n_components, component_mapping.size());
- const unsigned int n_threads = multithread_info.n_threads();
- Threads::ThreadGroup<> threads;
-
- // define starting and end point
- // for each thread
- typedef typename DoFHandler<dim,spacedim>::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::Mutex mutex;
-
- typedef std_cxx1x::tuple<const Mapping<dim,spacedim> &,
- const DoFHandler<dim,spacedim> &,
- const Quadrature<dim-1>&> Commons;
-
- // then assemble in parallel
- typedef void (*create_boundary_mass_matrix_1_t)
- (Commons,
- 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,
- const std::vector<unsigned int> &component_mapping,
- const MatrixCreator::internal::IteratorRange<DoFHandler<dim,spacedim> > range,
- Threads::Mutex &mutex);
- create_boundary_mass_matrix_1_t p
- = &create_boundary_mass_matrix_1<dim,spacedim>;
-
-//TODO: Use WorkStream here
- for (unsigned int thread=0; thread<n_threads; ++thread)
- threads += Threads::new_thread (p,
- Commons(mapping, dof, q), matrix,
- boundary_functions, rhs_vector,
- dof_to_boundary_mapping, coefficient,
- component_mapping,
- thread_ranges[thread], mutex);
- threads.join_all ();
+ MatrixCreator::internal::AssemblerBoundary::Scratch scratch;
+ MatrixCreator::internal::AssemblerBoundary::CopyData<dim,spacedim> copy_data;
+
+ WorkStream::run(dof.begin_active(),dof.end(),
+ static_cast<std_cxx1x::function<void (typename DoFHandler<dim,spacedim>::active_cell_iterator
+ const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<dim,spacedim> &)> >
+ (std_cxx1x::bind(create_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,std_cxx1x::_2,
+ std_cxx1x::_3,
+ std_cxx1x::cref(mapping),std_cxx1x::cref(fe),std_cxx1x::cref(q),
+ std_cxx1x::cref(boundary_functions),coefficient,
+ std_cxx1x::cref(component_mapping))),
+ static_cast<std_cxx1x::function<void (MatrixCreator::internal::AssemblerBoundary
+ ::CopyData<dim,spacedim> const &)> > (std_cxx1x::bind(
+ copy_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,
+ std_cxx1x::cref(boundary_functions),std_cxx1x::cref(dof_to_boundary_mapping),
+ std_cxx1x::ref(matrix),std_cxx1x::ref(rhs_vector))),
+ scratch,copy_data);
}