// <code>multithread_info</code> of that type) which can be used to query the
// number of processors in your system, which is often useful when deciding
// how many threads to start in parallel.
-#include <deal.II/base/thread_management.h>
+#include <deal.II/base/work_stream.h>
#include <deal.II/base/multithread_info.h>
// The next new include file declares a base class <code>TensorFunction</code>
{
using namespace dealii;
+ namespace Assembler
+ {
+ struct Scratch
+ {
+ Scratch() {};
+ };
+
+ struct CopyData
+ {
+ CopyData() {};
+
+ unsigned int dofs_per_cell;
+ std::vector<types::global_dof_index> local_dof_indices;
+ // We declare cell matrix and cell right hand side...
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+ };
+ }
+
// @sect3{AdvectionProblem class declaration}
// Following we declare the main class of this program. It is very much
// in the module, there are other, and possibly better suited, ways to
// achieve the same goal.
void assemble_system ();
- void assemble_system_interval (const typename DoFHandler<dim>::active_cell_iterator &begin,
- const typename DoFHandler<dim>::active_cell_iterator &end);
+ void build_local_system (typename DoFHandler<dim>::active_cell_iterator const &cell,
+ Assembler::Scratch &scratch,Assembler::CopyData ©_data);
+ void copy_local_to_global (Assembler::CopyData const ©_data);
+
// The following functions again are as in previous examples, as are the
// subsequent variables.
Vector<double> solution;
Vector<double> system_rhs;
-
- // When assembling the matrix in parallel, we have to synchronize when
- // several threads attempt to write the local contributions of a cell to
- // the global matrix at the same time. This is done using a
- // <code>Mutex</code>, which is an object that can be owned by only one
- // thread at a time. If a thread wants to write to the matrix, it has to
- // acquire this lock (if it is presently owned by another thread, then it
- // has to wait), then write to the matrix and finally release the
- // lock. Note that if the library was not compiled to support
- // multithreading (which you have to specify at the time you call the
- // <code>./configure</code> script in the top-level directory), then
- // the actual data type of the typedef
- // <code>Threads::Mutex</code> is a dummy class that provides all the
- // functions needed for a mutex, but does nothing when they are called;
- // this is reasonable, of course, since if only one thread is running at a
- // time, there is no need to synchronize with other threads.
- Threads::Mutex assembler_lock;
};
// is returned by n_threads(). This
// is also queried by functions inside the library to determine
// how many threads they shall create.
- const unsigned int n_threads = multithread_info.n_threads();
+
// It is worth noting, however, that this setup determines the load
// distribution onto processor in a static way: it does not take into
// account that some other part of our program may also be running
// object. Likewise, the function <code>join</code> that is supposed to
// wait for all spawned threads to return, returns immediately, as there
// can't be any threads running.
- Threads::ThreadGroup<> threads;
// Now we have to split the range of cells into chunks of approximately
// the same size. Each thread will then assemble the local contributions
// <code>raw_iterator</code>), and in this case the C++ language requires
// us to specify the template type explicitly. For brevity, we first
// typedef this data type to an alias.
+
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_handler.begin_active (),
- dof_handler.end (),
- n_threads);
// Finally, for each of the chunks of iterators we have computed, start
// one thread (or if not in multithread mode: execute assembly on these
// chunks sequentially). This is done using the following sequence of
// function calls:
- for (unsigned int thread=0; thread<n_threads; ++thread)
- threads += Threads::new_thread (&AdvectionProblem<dim>::assemble_system_interval,
- *this,
- thread_ranges[thread].first,
- thread_ranges[thread].second);
+
+ Assembler::Scratch scratch;
+ Assembler::CopyData copy_data;
+ WorkStream::run(dof_handler.begin_active(),dof_handler.end(),*this,
+ &AdvectionProblem::build_local_system,&AdvectionProblem::copy_local_to_global,
+ scratch,copy_data);
+
+
// The reasons and internal workings of these functions can be found in
// the report on the subject of multithreading, which is available online
// as well. Suffice it to say that we create a new thread that calls the
//
// Again, if the library was not configured to use multithreading, then
// no threads can run in parallel and the function returns immediately.
- threads.join_all ();
// After the matrix has been assembled in parallel, we still have to
template <int dim>
void
AdvectionProblem<dim>::
- assemble_system_interval (const typename DoFHandler<dim>::active_cell_iterator &begin,
- const typename DoFHandler<dim>::active_cell_iterator &end)
+ build_local_system (typename DoFHandler<dim>::active_cell_iterator const &cell,
+ Assembler::Scratch &scratch,
+ Assembler::CopyData ©_data)
{
// First of all, we will need some objects that describe boundary values,
// right hand side function and the advection field. As we will only
update_JxW_values | update_normal_vectors);
// Then we define some abbreviations to avoid unnecessarily long lines:
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ copy_data.dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
// We declare cell matrix and cell right hand side...
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_rhs (dofs_per_cell);
+ copy_data.cell_matrix = FullMatrix<double> (copy_data.dofs_per_cell, copy_data.dofs_per_cell);
+ copy_data.cell_rhs = Vector<double> (copy_data.dofs_per_cell);
// ... an array to hold the global indices of the degrees of freedom of
// the cell on which we are presently working...
- std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+ copy_data.local_dof_indices.resize(copy_data.dofs_per_cell);
// ... and array in which the values of right hand side, advection
// direction, and boundary values will be stored, for cell and face
std::vector<double> face_boundary_values (n_face_q_points);
std::vector<Tensor<1,dim> > face_advection_directions (n_face_q_points);
- // Then we start the main loop over the cells:
- typename DoFHandler<dim>::active_cell_iterator cell;
- for (cell=begin; cell!=end; ++cell)
- {
- // First clear old contents of the cell contributions...
- cell_matrix = 0;
- cell_rhs = 0;
-
- // ... then initialize the <code>FEValues</code> object...
- fe_values.reinit (cell);
-
- // ... obtain the values of right hand side and advection directions
- // at the quadrature points...
- advection_field.value_list (fe_values.get_quadrature_points(),
- advection_directions);
- right_hand_side.value_list (fe_values.get_quadrature_points(),
- rhs_values);
-
- // ... set the value of the streamline diffusion parameter as
- // described in the introduction...
- const double delta = 0.1 * cell->diameter ();
-
- // ... and assemble the local contributions to the system matrix and
- // right hand side as also discussed above:
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += ((advection_directions[q_point] *
- fe_values.shape_grad(j,q_point) *
- (fe_values.shape_value(i,q_point) +
- delta *
- (advection_directions[q_point] *
- fe_values.shape_grad(i,q_point)))) *
- fe_values.JxW(q_point));
-
- cell_rhs(i) += ((fe_values.shape_value(i,q_point) +
- delta *
- (advection_directions[q_point] *
- fe_values.shape_grad(i,q_point)) ) *
- rhs_values[q_point] *
- fe_values.JxW (q_point));
- };
- // Besides the cell terms which we have build up now, the bilinear
- // form of the present problem also contains terms on the boundary of
- // the domain. Therefore, we have to check whether any of the faces of
- // this cell are on the boundary of the domain, and if so assemble the
- // contributions of this face as well. Of course, the bilinear form
- // only contains contributions from the <code>inflow</code> part of
- // the boundary, but to find out whether a certain part of a face of
- // the present cell is part of the inflow boundary, we have to have
- // information on the exact location of the quadrature points and on
- // the direction of flow at this point; we obtain this information
- // using the FEFaceValues object and only decide within the main loop
- // whether a quadrature point is on the inflow boundary.
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face(face)->at_boundary())
- {
- // Ok, this face of the present cell is on the boundary of the
- // domain. Just as for the usual FEValues object which we have
- // used in previous examples and also above, we have to
- // reinitialize the FEFaceValues object for the present face:
- fe_face_values.reinit (cell, face);
-
- // For the quadrature points at hand, we ask for the values of
- // the inflow function and for the direction of flow:
- boundary_values.value_list (fe_face_values.get_quadrature_points(),
- face_boundary_values);
- advection_field.value_list (fe_face_values.get_quadrature_points(),
- face_advection_directions);
-
- // Now loop over all quadrature points and see whether it is on
- // the inflow or outflow part of the boundary. This is
- // determined by a test whether the advection direction points
- // inwards or outwards of the domain (note that the normal
- // vector points outwards of the cell, and since the cell is at
- // the boundary, the normal vector points outward of the domain,
- // so if the advection direction points into the domain, its
- // scalar product with the normal vector must be negative):
- for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
- if (fe_face_values.normal_vector(q_point) *
- face_advection_directions[q_point]
- < 0)
- // If the is part of the inflow boundary, then compute the
- // contributions of this face to the global matrix and right
- // hand side, using the values obtained from the
- // FEFaceValues object and the formulae discussed in the
- // introduction:
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) -= (face_advection_directions[q_point] *
- fe_face_values.normal_vector(q_point) *
- fe_face_values.shape_value(i,q_point) *
- fe_face_values.shape_value(j,q_point) *
- fe_face_values.JxW(q_point));
-
- cell_rhs(i) -= (face_advection_directions[q_point] *
- fe_face_values.normal_vector(q_point) *
- face_boundary_values[q_point] *
- fe_face_values.shape_value(i,q_point) *
- fe_face_values.JxW(q_point));
- };
- };
+ // ... then initialize the <code>FEValues</code> object...
+ fe_values.reinit (cell);
+
+ // ... obtain the values of right hand side and advection directions
+ // at the quadrature points...
+ advection_field.value_list (fe_values.get_quadrature_points(),
+ advection_directions);
+ right_hand_side.value_list (fe_values.get_quadrature_points(),
+ rhs_values);
+
+ // ... set the value of the streamline diffusion parameter as
+ // described in the introduction...
+ const double delta = 0.1 * cell->diameter ();
+
+ // ... and assemble the local contributions to the system matrix and
+ // right hand side as also discussed above:
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+ copy_data.cell_matrix(i,j) += ((advection_directions[q_point] *
+ fe_values.shape_grad(j,q_point) *
+ (fe_values.shape_value(i,q_point) +
+ delta *
+ (advection_directions[q_point] *
+ fe_values.shape_grad(i,q_point)))) *
+ fe_values.JxW(q_point));
+
+ copy_data.cell_rhs(i) += ((fe_values.shape_value(i,q_point) +
+ delta *
+ (advection_directions[q_point] *
+ fe_values.shape_grad(i,q_point)) ) *
+ rhs_values[q_point] *
+ fe_values.JxW (q_point));
+ };
+
+ // Besides the cell terms which we have build up now, the bilinear
+ // form of the present problem also contains terms on the boundary of
+ // the domain. Therefore, we have to check whether any of the faces of
+ // this cell are on the boundary of the domain, and if so assemble the
+ // contributions of this face as well. Of course, the bilinear form
+ // only contains contributions from the <code>inflow</code> part of
+ // the boundary, but to find out whether a certain part of a face of
+ // the present cell is part of the inflow boundary, we have to have
+ // information on the exact location of the quadrature points and on
+ // the direction of flow at this point; we obtain this information
+ // using the FEFaceValues object and only decide within the main loop
+ // whether a quadrature point is on the inflow boundary.
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->at_boundary())
+ {
+ // Ok, this face of the present cell is on the boundary of the
+ // domain. Just as for the usual FEValues object which we have
+ // used in previous examples and also above, we have to
+ // reinitialize the FEFaceValues object for the present face:
+ fe_face_values.reinit (cell, face);
+
+ // For the quadrature points at hand, we ask for the values of
+ // the inflow function and for the direction of flow:
+ boundary_values.value_list (fe_face_values.get_quadrature_points(),
+ face_boundary_values);
+ advection_field.value_list (fe_face_values.get_quadrature_points(),
+ face_advection_directions);
+
+ // Now loop over all quadrature points and see whether it is on
+ // the inflow or outflow part of the boundary. This is
+ // determined by a test whether the advection direction points
+ // inwards or outwards of the domain (note that the normal
+ // vector points outwards of the cell, and since the cell is at
+ // the boundary, the normal vector points outward of the domain,
+ // so if the advection direction points into the domain, its
+ // scalar product with the normal vector must be negative):
+ for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+ if (fe_face_values.normal_vector(q_point) *
+ face_advection_directions[q_point]
+ < 0)
+ // If the is part of the inflow boundary, then compute the
+ // contributions of this face to the global matrix and right
+ // hand side, using the values obtained from the
+ // FEFaceValues object and the formulae discussed in the
+ // introduction:
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+ copy_data.cell_matrix(i,j) -= (face_advection_directions[q_point] *
+ fe_face_values.normal_vector(q_point) *
+ fe_face_values.shape_value(i,q_point) *
+ fe_face_values.shape_value(j,q_point) *
+ fe_face_values.JxW(q_point));
+
+ copy_data.cell_rhs(i) -= (face_advection_directions[q_point] *
+ fe_face_values.normal_vector(q_point) *
+ face_boundary_values[q_point] *
+ fe_face_values.shape_value(i,q_point) *
+ fe_face_values.JxW(q_point));
+ };
+ };
+
+
+ // Now go on by transferring the local contributions to the system of
+ // equations into the global objects. The first step was to obtain the
+ // global indices of the degrees of freedom on this cell.
+ cell->get_dof_indices (copy_data.local_dof_indices);
+ }
- // Now go on by transferring the local contributions to the system of
- // equations into the global objects. The first step was to obtain the
- // global indices of the degrees of freedom on this cell.
- cell->get_dof_indices (local_dof_indices);
+
+ template <int dim>
+ void
+ AdvectionProblem<dim>::copy_local_to_global (Assembler::CopyData const ©_data)
+ {
// Up until now we have not taken care of the fact that this function
// might run more than once in parallel, as the operations above only
// matrix, and can do so freely. When finished, we release the lock
// again so as to allow other threads to acquire it and write to the
// matrix.
- assembler_lock.acquire ();
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
{
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- cell_matrix(i,j));
+ for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+ system_matrix.add (copy_data.local_dof_indices[i],
+ copy_data.local_dof_indices[j],
+ copy_data.cell_matrix(i,j));
- system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ system_rhs(copy_data.local_dof_indices[i]) += copy_data.cell_rhs(i);
};
- assembler_lock.release ();
// At this point, the locked operations on the global matrix are done,
// i.e. other threads can now enter into the protected section by
// acquiring the lock. Two final notes are in place here, however:
// class; since all threads execute member functions of the same
// object, they have the same <code>this</code> pointer and therefore
// also operate on the same <code>lock</code>.
- };
}
+
#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 <typename DH>
+ 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 DH::active_cell_iterator cell;
+ std::vector<FullMatrix<double> > cell_matrix;
+ std::vector<Vector<double> > cell_vector;
+ };
+
+ template <typename DH>
+ CopyData<DH>::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<DoFHandler<dim,
+ spacedim> > ©_data,
+ Mapping<dim, spacedim> const &mapping,
+ FiniteElement<dim,spacedim> const &fe,
+ Quadrature<dim-1> const &q,
+ typename FunctionMap<spacedim>::type const &boundary_functions,
+ Function<spacedim> const *const coefficient,
+ std::vector<unsigned int> const &component_mapping)
+
{
// All assertions for this function
// are in the calling function
// before creating threads.
- const 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 |
update_JxW_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));
-
- std::vector<types::global_dof_index> dofs (dofs_per_cell);
- std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face);
+ Vector<double>(n_function_components));
- // 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,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;
+ copy_data.dofs.resize(copy_data.dofs_per_cell);
+ cell->get_dof_indices (copy_data.dofs);
- fe_values.reinit (cell, face);
+ std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_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_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
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- coefficient_vector_values[point] = coefficient_values[point];
- }
+ // 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);
- // 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);
+
+ 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];
+ }
- 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;
+ // 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 (!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 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;
- 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);
+ 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);
+ }
- 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());
-
- // 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)
+ 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);
+
+ 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));
+ // check for each of the dofs on this cell
+ // whether it is on the face
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(),
+ dofs_on_face_vector.end(),
+ copy_data.dofs[i])
+ !=
+ dofs_on_face_vector.end());
+ }
+ }
+
+ template <int dim,int spacedim>
+ void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,
+ spacedim> > const ©_data,
+ typename FunctionMap<spacedim>::type const &boundary_functions,
+ std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
+ SparseMatrix<double> &matrix,
+ Vector<double> &rhs_vector)
+ {
+ // now transfer cell matrix and vector to the whole boundary matrix
+ //
+ // in the following: dof[i] holds the global index of the i-th degree of
+ // freedom on the present cell. If it is also a dof on the boundary, it
+ // must be a nonzero entry in the dof_to_boundary_mapping and then
+ // the boundary index of this dof is dof_to_boundary_mapping[dof[i]].
+ //
+ // if dof[i] is not on the boundary, it should be zero on the boundary
+ // therefore on all quadrature points and finally all of its
+ // entries in the cell matrix and vector should be zero. If not, we
+ // throw an error (note: because of the evaluation of the shape
+ // functions only up to machine precision, the term "must be zero"
+ // really should mean: "should be very small". since this is only an
+ // assertion and not part of the code, we may choose "very small"
+ // quite arbitrarily)
+ //
+ // the main problem here is that the matrix or vector entry should also
+ // be zero if the degree of freedom dof[i] is on the boundary, but not
+ // on the present face, i.e. on another face of the same cell also
+ // on the boundary. We can therefore not rely on the
+ // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to
+ // determine whether dof[i] is a dof on the present face. We do so
+ // by getting the dofs on the face into @p{dofs_on_face_vector},
+ // a vector as always. Usually, searching in a vector is
+ // inefficient, so we copy the dofs into a set, which enables binary
+ // searches.
+ unsigned int pos(0);
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ {
+ // check if this face is on that part of
+ // the boundary we are interested in
+ if (boundary_functions.find(copy_data.cell->face(face)->boundary_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> (DoFHandler<2,3>::active_cell_iterator const &cell,
+ MatrixCreator::internal::AssemblerBoundary::Scratch const
+ &scratch,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<2,
+ 3> > ©_data,
+ Mapping<2,3> const &mapping,
+ FiniteElement<2,3> const &fe,
+ Quadrature<1> const &q,
+ 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> (DoFHandler<1,3>::active_cell_iterator const &cell,
+ MatrixCreator::internal::AssemblerBoundary::Scratch const
+ &scratch,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<1,
+ 3> > ©_data,
+ Mapping<1,3> const &mapping,
+ FiniteElement<1,3> const &fe,
+ Quadrature<0> const &q,
+ 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<DoFHandler<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<DoFHandler<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<DoFHandler<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);
}
template <int dim, int spacedim>
void
- create_boundary_mass_matrix_1 (std_cxx1x::tuple<const hp::MappingCollection<dim,spacedim> &,
- const hp::DoFHandler<dim,spacedim> &,
- const hp::QCollection<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<hp::DoFHandler<dim,spacedim> > range,
- Threads::Mutex &mutex)
+ create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator const
+ &cell,
+ MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch,
+ MatrixCreator::internal::AssemblerBoundary
+ ::CopyData<hp::DoFHandler<dim,spacedim> > ©_data,
+ hp::MappingCollection<dim,spacedim> const &mapping,
+ hp::FECollection<dim,spacedim> const &fe_collection,
+ hp::QCollection<dim-1> const &q,
+ const typename FunctionMap<spacedim>::type &boundary_functions,
+ Function<spacedim> const *const coefficient,
+ std::vector<unsigned int> const &component_mapping)
{
- const hp::MappingCollection<dim,spacedim> &mapping = std_cxx1x::get<0>(commons);
- const hp::DoFHandler<dim,spacedim> &dof = std_cxx1x::get<1>(commons);
- const hp::QCollection<dim-1>& q = std_cxx1x::get<2>(commons);
- const hp::FECollection<dim,spacedim> &fe_collection = dof.get_fe();
const unsigned int n_components = fe_collection.n_components();
const unsigned int n_function_components = boundary_functions.begin()->second->n_components;
const bool fe_is_system = (n_components != 1);
-#ifdef DEBUG
- if (true)
- {
- types::global_dof_index max_element = static_cast<types::global_dof_index>(0);
- for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
- i!=dof_to_boundary_mapping.end(); ++i)
- if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) &&
- (*i > max_element))
- max_element = *i;
- Assert (max_element == matrix.n()-1, ExcInternalError());
- };
-#endif
+ const FiniteElement<dim,spacedim> &fe = cell->get_fe();
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+
+ copy_data.cell = cell;
+ copy_data.dofs_per_cell = fe.dofs_per_cell;
+ copy_data.dofs.resize(copy_data.dofs_per_cell);
+ cell->get_dof_indices (copy_data.dofs);
const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell(),
max_dofs_per_face = fe_collection.max_dofs_per_face();
std::vector<double> rhs_values_scalar;
std::vector<Vector<double> > rhs_values_system;
- std::vector<types::global_dof_index> dofs (max_dofs_per_cell);
- std::vector<types::global_dof_index> dofs_on_face_vector (max_dofs_per_face);
+ 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(max_dofs_per_cell);
+ copy_data.dofs.resize(copy_data.dofs_per_cell);
+ cell->get_dof_indices (copy_data.dofs);
+ // Because CopyData objects are reused and that push_back is used,
+ // dof_is_on_face, cell_matrix, and cell_vector must be cleared before
+ // they are reused
+ copy_data.dof_is_on_face.clear();
+ copy_data.cell_matrix.clear();
+ copy_data.cell_vector.clear();
- typename hp::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())
- {
- x_fe_values.reinit (cell, face);
- const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values ();
+ 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())
+ {
+ x_fe_values.reinit (cell, face);
- const FiniteElement<dim,spacedim> &fe = cell->get_fe();
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int dofs_per_face = fe.dofs_per_face;
+ const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values ();
- cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
- cell_vector.reinit (dofs_per_cell);
- cell_matrix = 0;
- cell_vector = 0;
+ copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
+ copy_data.dofs_per_cell));
+ copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
- if (fe_is_system)
- // FE has several components
- {
- rhs_values_system.resize (fe_values.n_quadrature_points,
- Vector<double>(n_function_components));
- boundary_functions.find(cell->face(face)->boundary_indicator())
- ->second->vector_value_list (fe_values.get_quadrature_points(),
- rhs_values_system);
+ if (fe_is_system)
+ // FE has several components
+ {
+ rhs_values_system.resize (fe_values.n_quadrature_points,
+ Vector<double>(n_function_components));
+ boundary_functions.find(cell->face(face)->boundary_indicator())
+ ->second->vector_value_list (fe_values.get_quadrature_points(),
+ rhs_values_system);
- if (coefficient != 0)
- {
- if (coefficient->n_components==1)
- {
- coefficient_values.resize (fe_values.n_quadrature_points);
- coefficient->value_list (fe_values.get_quadrature_points(),
- coefficient_values);
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
- {
- const double v = fe_values.shape_value(i,point);
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- if (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first)
- {
- const double u = fe_values.shape_value(j,point);
- cell_matrix(i,j)
- += (u * v * weight * coefficient_values[point]);
- }
-
- cell_vector(i) += v *
- rhs_values_system[point](
- component_mapping[fe.system_to_component_index(i).first]) * weight;
- }
- }
- }
- else
- {
- coefficient_vector_values.resize (fe_values.n_quadrature_points,
- Vector<double>(n_components));
- coefficient->vector_value_list (fe_values.get_quadrature_points(),
- coefficient_vector_values);
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
- {
- const double v = fe_values.shape_value(i,point);
- const unsigned int component_i=
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- if (fe.system_to_component_index(j).first ==
- component_i)
- {
- const double u = fe_values.shape_value(j,point);
- cell_matrix(i,j) +=
- (u * v * weight * coefficient_vector_values[point](component_i));
- }
- cell_vector(i) += v * rhs_values_system[point](component_mapping[component_i]) * weight;
- }
- }
- }
- }
- else // if (coefficient == 0)
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ if (coefficient != 0)
+ {
+ if (coefficient->n_components==1)
{
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+ coefficient_values.resize (fe_values.n_quadrature_points);
+ coefficient->value_list (fe_values.get_quadrature_points(),
+ coefficient_values);
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
{
- const double v = fe_values.shape_value(i,point);
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- if (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first)
- {
- const double u = fe_values.shape_value(j,point);
- cell_matrix(i,j) += (u * v * weight);
- }
- cell_vector(i) += v *
- rhs_values_system[point](
- fe.system_to_component_index(i).first) *
- weight;
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+ {
+ const double v = fe_values.shape_value(i,point);
+ for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+ if (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first)
+ {
+ const double u = fe_values.shape_value(j,point);
+ copy_data.cell_matrix.back()(i,j)
+ += (u * v * weight * coefficient_values[point]);
+ }
+
+ copy_data.cell_vector.back()(i) += v *
+ rhs_values_system[point](
+ component_mapping[fe.system_to_component_index(i).first]) * weight;
+ }
}
}
- }
- else
- // FE is a scalar one
- {
- rhs_values_scalar.resize (fe_values.n_quadrature_points);
- boundary_functions.find(cell->face(face)->boundary_indicator())
- ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
-
- if (coefficient != 0)
+ else
+ {
+ coefficient_vector_values.resize (fe_values.n_quadrature_points,
+ Vector<double>(n_components));
+ coefficient->vector_value_list (fe_values.get_quadrature_points(),
+ coefficient_vector_values);
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+ {
+ const double v = fe_values.shape_value(i,point);
+ const unsigned int component_i=
+ fe.system_to_component_index(i).first;
+ for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+ if (fe.system_to_component_index(j).first ==
+ component_i)
+ {
+ const double u = fe_values.shape_value(j,point);
+ copy_data.cell_matrix.back()(i,j) +=
+ (u * v * weight * coefficient_vector_values[point](component_i));
+ }
+ copy_data.cell_vector.back()(i) += v *
+ rhs_values_system[point](component_mapping[component_i]) * weight;
+ }
+ }
+ }
+ }
+ else // if (coefficient == 0)
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
{
- coefficient_values.resize (fe_values.n_quadrature_points);
- coefficient->value_list (fe_values.get_quadrature_points(),
- coefficient_values);
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
{
- const double 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;
- }
+ const double v = fe_values.shape_value(i,point);
+ for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+ if (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first)
+ {
+ const double u = fe_values.shape_value(j,point);
+ copy_data.cell_matrix.back()(i,j) += (u * v * weight);
+ }
+ copy_data.cell_vector.back()(i) += v *
+ rhs_values_system[point](
+ fe.system_to_component_index(i).first) *
+ weight;
}
}
- else
+ }
+ else
+ // FE is a scalar one
+ {
+ rhs_values_scalar.resize (fe_values.n_quadrature_points);
+ boundary_functions.find(cell->face(face)->boundary_indicator())
+ ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
+
+ if (coefficient != 0)
+ {
+ coefficient_values.resize (fe_values.n_quadrature_points);
+ coefficient->value_list (fe_values.get_quadrature_points(),
+ coefficient_values);
for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
{
const double weight = fe_values.JxW(point);
for (unsigned int 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);
+ copy_data.cell_matrix.back()(i,j) += (u * v * weight *
+ coefficient_values[point]);
}
- cell_vector(i) += v * rhs_values_scalar[point] * weight;
+ copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight;
}
}
- }
+ }
+ else
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+ {
+ const double v = fe_values.shape_value(i,point);
+ for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+ {
+ const double u = fe_values.shape_value(j,point);
+ copy_data.cell_matrix.back()(i,j) += (u * v * weight);
+ }
+ copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] * weight;
+ }
+ }
+ }
+
+ cell->face(face)->get_dof_indices (dofs_on_face_vector,
+ cell->active_fe_index());
+ // for each dof on the cell, have a
+ // flag whether it is on the face
+ copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell));
+ // check for each of the dofs on this cell
+ // whether it is on the face
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(),
+ dofs_on_face_vector.end(),
+ copy_data.dofs[i])
+ !=
+ dofs_on_face_vector.end());
+ }
+ }
+
- // 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.
- dofs.resize (dofs_per_cell);
- dofs_on_face_vector.resize (dofs_per_face);
- dof_is_on_face.resize (dofs_per_cell);
-
- cell->get_dof_indices (dofs);
- cell->face(face)->get_dof_indices (dofs_on_face_vector,
- cell->active_fe_index());
-
- // 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.
+
+ template <int dim,int spacedim>
+ void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary
+ ::CopyData<hp::DoFHandler<dim,spacedim> > const ©_data,
+ typename FunctionMap<spacedim>::type const &boundary_functions,
+ std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
+ SparseMatrix<double> &matrix,
+ Vector<double> &rhs_vector)
+ {
+ // now transfer cell matrix and vector to the whole boundary matrix
+ //
+ // in the following: dof[i] holds the global index of the i-th degree of
+ // freedom on the present cell. If it is also a dof on the boundary, it
+ // must be a nonzero entry in the dof_to_boundary_mapping and then
+ // the boundary index of this dof is dof_to_boundary_mapping[dof[i]].
+ //
+ // if dof[i] is not on the boundary, it should be zero on the boundary
+ // therefore on all quadrature points and finally all of its
+ // entries in the cell matrix and vector should be zero. If not, we
+ // throw an error (note: because of the evaluation of the shape
+ // functions only up to machine precision, the term "must be zero"
+ // really should mean: "should be very small". since this is only an
+ // assertion and not part of the code, we may choose "very small"
+ // quite arbitrarily)
+ //
+ // the main problem here is that the matrix or vector entry should also
+ // be zero if the degree of freedom dof[i] is on the boundary, but not
+ // on the present face, i.e. on another face of the same cell also
+ // on the boundary. We can therefore not rely on the
+ // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to
+ // determine whether dof[i] is a dof on the present face. We do so
+ // by getting the dofs on the face into @p{dofs_on_face_vector},
+ // a vector as always. Usually, searching in a vector is
+ // inefficient, so we copy the dofs into a set, which enables binary
+ // searches.
+ unsigned int pos(0);
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ {
+ // check if this face is on that part of
+ // the boundary we are interested in
+ if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) !=
+ boundary_functions.end())
+ {
#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));
+ // in debug mode: compute an element in the matrix which is
+ // guaranteed to belong to a boundary dof. We do this to check that the
+ // entries in the cell matrix are guaranteed to be zero if the
+ // respective dof is not on the boundary. Since because of
+ // round-off, the actual value of the matrix entry may be
+ // only close to zero, we assert that it is small relative to an element
+ // which is guaranteed to be nonzero. (absolute smallness does not
+ // suffice since the size of the domain scales in here)
+ //
+ // for this purpose we seek the diagonal of the matrix, where there
+ // must be an element belonging to the boundary. we take the maximum
+ // diagonal entry.
+ types::global_dof_index max_element = static_cast<types::global_dof_index>(0);
+ for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
+ i!=dof_to_boundary_mapping.end(); ++i)
+ if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) &&
+ (*i > max_element))
+ max_element = *i;
+ Assert (max_element == matrix.n()-1, ExcInternalError());
+
+ double max_diag_entry = 0;
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ if (std::fabs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry)
+ max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i));
#endif
- // lock the matrix
- Threads::Mutex::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);
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+ {
+ if (copy_data.dof_is_on_face[pos][i] && copy_data.dof_is_on_face[pos][j])
+ matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
+ dof_to_boundary_mapping[copy_data.dofs[j]],
+ copy_data.cell_matrix[pos](i,j));
else
{
- // compare here for relative
- // smallness
- Assert (std::fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
- ExcInternalError());
+ // assume that all shape functions that are nonzero on the boundary
+ // are also listed in the @p{dof_to_boundary} mapping. if that
+ // is not the case, then the boundary mass matrix does not
+ // make that much sense anyway, as it only contains entries for
+ // parts of the functions living on the boundary
+ //
+ // these, we may compare here for relative smallness of all
+ // entries in the local matrix which are not taken over to
+ // the global one
+ Assert (std::fabs(copy_data.cell_matrix[pos](i,j)) <= 1e-10 * max_diag_entry,
+ ExcInternalError ());
}
- }
- }
+ }
+ for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+ if (copy_data.dof_is_on_face[pos][j])
+ rhs_vector(dof_to_boundary_mapping[copy_data.dofs[j]]) += copy_data.cell_vector[pos](j);
+ else
+ {
+ // compare here for relative
+ // smallness
+ Assert (std::fabs(copy_data.cell_vector[pos](j)) <= 1e-10 * max_diag_entry,
+ ExcInternalError());
+ }
+ ++pos;
+ }
+ }
+ }
}
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 hp::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);
-
- typedef std_cxx1x::tuple<const hp::MappingCollection<dim,spacedim> &,
- const hp::DoFHandler<dim,spacedim> &,
- const hp::QCollection<dim-1>&> Commons;
-
- // mutex to synchronise access to
- // the matrix
- Threads::Mutex mutex;
-
- // 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<hp::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<hp::DoFHandler<dim,spacedim> > copy_data;
+
+ WorkStream::run(dof.begin_active(),dof.end(),
+ static_cast<std_cxx1x::function<void (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
+ const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > &)> >
+ (std_cxx1x::bind(create_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,std_cxx1x::_2,
+ std_cxx1x::_3,
+ std_cxx1x::cref(mapping),std_cxx1x::cref(fe_collection),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<hp::DoFHandler<dim,spacedim> > const &)> > (std_cxx1x::bind(
+ copy_hp_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);
}
projects / dealii-svn.git / commitdiff