From: heltai Date: Thu, 17 Oct 2013 15:29:01 +0000 (+0000) Subject: Merged from trunk X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f2b21ce0d3791c91465944f0a5d4e35dce6df743;p=dealii-svn.git Merged from trunk git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@31284 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/cmake/setup_compiler_flags_intel.cmake b/deal.II/cmake/setup_compiler_flags_intel.cmake index 912798af47..98afb1943b 100644 --- a/deal.II/cmake/setup_compiler_flags_intel.cmake +++ b/deal.II/cmake/setup_compiler_flags_intel.cmake @@ -62,6 +62,14 @@ ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-w2") # -w68 integer conversion resulted in a change of sign # (triggers a lot in functionparser) # -w175 subscript out of range +# -w135 class template "dealii::FE_Q_Base" +# has no member "Implementation" +# (the compiler is objectively wrong since the warning +# triggers also on code of the form +# class FE_Q_Base { +# struct Implementation; // forward declaration +# friend struct Implementation; +# };) # -w177 declared but not referenced # -w279 controlling expression is constant # -w327 NULL reference is not allowed @@ -78,6 +86,7 @@ ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-w2") # -w2536 type qualifiers are meaningless here # ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd68") +ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd135") ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd175") ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd177") ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd279") diff --git a/deal.II/doc/doxygen/headers/global_dof_index.h b/deal.II/doc/doxygen/headers/global_dof_index.h index 51ebce36d7..2b5d84bddd 100644 --- a/deal.II/doc/doxygen/headers/global_dof_index.h +++ b/deal.II/doc/doxygen/headers/global_dof_index.h @@ -18,10 +18,22 @@ /** * @page GlobalDoFIndex When to use types::global_dof_index instead of unsigned int * - * When the 64-bits version of deal.II is used, it becomes necessary to - * declare as types::global_dof_index (unsigned long long int) instead of - * unsigned int. Here, we want to clarify when types::global_dof_index must be - * used. + * deal.II can be configured to use 64-bit indices for degrees of freedom, + * rather than the usual unsigned integers that default to 32-bit on most + * current systems. This is necessary since we want to be able to solve + * problems with more than 4 billion unknowns (the limit of what can be + * represented with 32-bit unsigned integers). At the same time, we do not + * want to indiscriminately replace all integers in deal.II with 64-bit + * versions, since this would increase memory use in many places where we + * represent quantities that will most definitely not be larger than 4 billion. + * + * The data type we define for these indices to keep the bulk + * of the code base free of #ifdefs is types::global_dof_index. + * If deal.II is configured as normal, this type is unsigned int, + * but can be switched to unsigned long long int if the right + * flag is provided (see the ReadMe file). This page is intended to clarify + * when types::global_dof_index must be used and when one can use a regular + * unsigned integer: * *
* @@ -35,15 +47,18 @@ * *
@anchor GlobalDoFIndexCell Cell
*
- * The ID of cell is not unique. Cells with different levels of refinement + * The ID of cell is not unique: Cells with different levels of refinement * and/or on different processors can have the same ID. Thus, all the data - * associated to cells can be unsigned int. + * associated to cells can be unsigned int because on a single processor, + * one one mesh level, there will definitely not be more than 4 billion + * cells. *
* *
@anchor GlobalDoFIndexDoFHandler * DoFHandler
*
- * The ID of each degree of freedom is unique. Therefore, degrees of freedom + * The ID of each degree of freedom is unique in a parallel computation. + * Therefore, degrees of freedom * are types::global_dof_index. *
* @@ -52,7 +67,8 @@ *
* The numbers of row and column are types::global_dof_index even if it is not * expected that someone will create a FullMatrix with so many entries. - * However, ConstraintMatrix is templated on the matrix type and thus, the + * However, some functions of ConstraintMatrix are templated on the matrix + * type and thus, the * size of a FullMatrix has to be of the same type than the size of * SparseMatrix. *
@@ -60,11 +76,13 @@ *
@anchor GlobalDoFIndexSparseMatrix * SparseMatrix
*
- * The size of SparseMatrix can be arbitrary large therefore, - * types::global_do_index is used. However, even for a large complex problem we - * can solve now, there is no reason for the number of non-zero entries in a - * sparse matrix to go over four billions. Thus, we still use unsigned int - * for, e.g., row_lengths in the object. + * The size of SparseMatrix can be arbitrary large and it is conceivable that + * with sufficient memory on a single node, one may generate a matrix with + * more than 4 billion rows or columns. Therefore, types::global_dof_index is + * used. However, even for a large complex problem we can solve now, it is not + * reasonable to expect the number of non-zero entries in a sparse matrix to + * go over four billion. Thus, we still use unsigned int for, e.g., + * SparsityPattern::row_lengths and similar functions. *
* *
diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 7cdcdd827e..54d0a2412b 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -145,6 +145,12 @@ inconvenience this causes.

Specific improvements

    +
  1. + Fixed: When deriving from DataOut to filter the cells where output is generated, there were two different bugs that result in segmentation faults or wrong cells written (example, step-18). +
    + (Timo Heister, 2013/10/16) +
    2. +
  2. New: GridIn::read_vtk() reads 2d and 3d meshes in VTK format.
    diff --git a/deal.II/examples/step-20/doc/intro.dox b/deal.II/examples/step-20/doc/intro.dox index dd4c3829ca..13476bc201 100644 --- a/deal.II/examples/step-20/doc/intro.dox +++ b/deal.II/examples/step-20/doc/intro.dox @@ -400,10 +400,10 @@ also a full matrix. On the other hand, the CG algorithm doesn't require us to actually have a representation of $S$, it is sufficient to form matrix-vector products with it. We can do so in steps: to compute $Sv=B^TM^{-1}Bv=B^T(M^{-1}(Bv))$, we
      -
    1. form $w = T v$; +
    2. form $w = B v$;
    3. solve $My = w$ for $y=M^{-1}w$, using the CG method applied to the positive definite and symmetric mass matrix $M$; -
    4. form $z=B^Ty$ to obtain $Sv=z$. +
    5. form $z=B^Ty$ to obtain $z=Sv$.
    Note how we evaluate the expression $B^TM^{-1}Bv$ right to left to avoid matrix-matrix products; this way, all we have to do is evaluate diff --git a/deal.II/examples/step-9/step-9.cc b/deal.II/examples/step-9/step-9.cc index 4a9e28152e..73b56ce912 100644 --- a/deal.II/examples/step-9/step-9.cc +++ b/deal.II/examples/step-9/step-9.cc @@ -54,7 +54,7 @@ // multithread_info 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 +#include #include // The next new include file declares a base class TensorFunction @@ -74,6 +74,25 @@ namespace Step9 { using namespace dealii; + namespace Assembler + { + struct Scratch + { + Scratch() {}; + }; + + struct CopyData + { + CopyData() {}; + + unsigned int dofs_per_cell; + std::vector local_dof_indices; + // We declare cell matrix and cell right hand side... + FullMatrix cell_matrix; + Vector cell_rhs; + }; + } + // @sect3{AdvectionProblem class declaration} // Following we declare the main class of this program. It is very much @@ -107,8 +126,10 @@ namespace Step9 // 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::active_cell_iterator &begin, - const typename DoFHandler::active_cell_iterator &end); + void build_local_system (typename DoFHandler::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. @@ -128,23 +149,6 @@ namespace Step9 Vector solution; Vector 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 - // Mutex, 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 - // ./configure script in the top-level directory), then - // the actual data type of the typedef - // Threads::Mutex 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; }; @@ -526,7 +530,7 @@ namespace Step9 // 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 @@ -559,7 +563,6 @@ namespace Step9 // object. Likewise, the function join 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 @@ -584,22 +587,21 @@ namespace Step9 // raw_iterator), 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::active_cell_iterator active_cell_iterator; - std::vector > - thread_ranges - = Threads::split_range (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::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 @@ -624,7 +626,6 @@ namespace Step9 // // 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 @@ -647,8 +648,9 @@ namespace Step9 template void AdvectionProblem:: - assemble_system_interval (const typename DoFHandler::active_cell_iterator &begin, - const typename DoFHandler::active_cell_iterator &end) + build_local_system (typename DoFHandler::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 @@ -681,17 +683,17 @@ namespace Step9 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 cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); + copy_data.cell_matrix = FullMatrix (copy_data.dofs_per_cell, copy_data.dofs_per_cell); + copy_data.cell_rhs = Vector (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 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 @@ -701,117 +703,118 @@ namespace Step9 std::vector face_boundary_values (n_face_q_points); std::vector > face_advection_directions (n_face_q_points); - // Then we start the main loop over the cells: - typename DoFHandler::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 FEValues 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_pointinflow 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::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_pointFEValues 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_pointinflow 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::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_pointget_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 + void + AdvectionProblem::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 @@ -836,17 +839,15 @@ namespace Step9 // 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; ithis pointer and therefore // also operate on the same lock. - }; } + diff --git a/deal.II/include/deal.II/base/types.h b/deal.II/include/deal.II/base/types.h index 0da2ba52e6..7fa251bcad 100644 --- a/deal.II/include/deal.II/base/types.h +++ b/deal.II/include/deal.II/base/types.h @@ -70,6 +70,9 @@ namespace types * * The data type always indicates an * unsigned integer type. + * + * See the @ref GlobalDoFIndex page for guidance on when this + * type should or should not be used. */ // TODO: we should check that unsigned long long int // has the same size as uint64_t diff --git a/deal.II/include/deal.II/base/vectorization.h b/deal.II/include/deal.II/base/vectorization.h index f2c560f36c..ec48938935 100644 --- a/deal.II/include/deal.II/base/vectorization.h +++ b/deal.II/include/deal.II/base/vectorization.h @@ -1952,10 +1952,8 @@ DEAL_II_NAMESPACE_CLOSE namespace std { /** - * Computes the sine of a vectorized - * data field. The result is return as - * vectorized array in the form - * {sin(x[0]), sin(x[1]), ..., + * Computes the sine of a vectorized data field. The result is return as + * vectorized array in the form {sin(x[0]), sin(x[1]), ..., * sin(x[n_array_elements-1])}. * * @relates VectorizedArray @@ -1974,10 +1972,8 @@ namespace std /** - * Computes the tangent of a vectorized - * data field. The result is return as - * vectorized array in the form - * {tan(x[0]), tan(x[1]), ..., + * Computes the tangent of a vectorized data field. The result is return as + * vectorized array in the form {tan(x[0]), tan(x[1]), ..., * tan(x[n_array_elements-1])}. * * @relates VectorizedArray @@ -1995,10 +1991,8 @@ namespace std /** - * Computes the cosine of a vectorized - * data field. The result is return as - * vectorized array in the form - * {cos(x[0]), cos(x[1]), ..., + * Computes the cosine of a vectorized data field. The result is return as + * vectorized array in the form {cos(x[0]), cos(x[1]), ..., * cos(x[n_array_elements-1])}. * * @relates VectorizedArray @@ -2016,10 +2010,8 @@ namespace std /** - * Computes the exponential of a vectorized - * data field. The result is return as - * vectorized array in the form - * {exp(x[0]), exp(x[1]), ..., + * Computes the exponential of a vectorized data field. The result is return + * as vectorized array in the form {exp(x[0]), exp(x[1]), ..., * exp(x[n_array_elements-1])}. * * @relates VectorizedArray @@ -2037,10 +2029,8 @@ namespace std /** - * Computes the natural logarithm of a - * vectorized data field. The result is return - * as vectorized array in the form - * {log(x[0]), log(x[1]), ..., + * Computes the natural logarithm of a vectorized data field. The result is + * return as vectorized array in the form {log(x[0]), log(x[1]), ..., * log(x[n_array_elements-1])}. * * @relates VectorizedArray @@ -2059,10 +2049,8 @@ namespace std /** - * Computes the square root of a vectorized - * data field. The result is return as - * vectorized array in the form - * {sqrt(x[0]), sqrt(x[1]), ..., + * Computes the square root of a vectorized data field. The result is return + * as vectorized array in the form {sqrt(x[0]), sqrt(x[1]), ..., * sqrt(x[n_array_elements-1])}. * * @relates VectorizedArray @@ -2078,11 +2066,9 @@ namespace std /** - * Computes the absolute value (modulus) of a - * vectorized data field. The result is return - * as vectorized array in the form - * {abs(x[0]), abs(x[1]), ..., - * abs(x[n_array_elements-1])}. + * Computes the absolute value (modulus) of a vectorized data field. The + * result is return as vectorized array in the form {abs(x[0]), + * abs(x[1]), ..., abs(x[n_array_elements-1])}. * * @relates VectorizedArray */ @@ -2097,11 +2083,9 @@ namespace std /** - * Computes the componentwise maximum of two - * vectorized data fields. The result is - * return as vectorized array in the form - * {max(x[0],y[0]), max(x[1],y[1]), - * ...}. + * Computes the componentwise maximum of two vectorized data fields. The + * result is return as vectorized array in the form {max(x[0],y[0]), + * max(x[1],y[1]), ...}. * * @relates VectorizedArray */ @@ -2117,11 +2101,9 @@ namespace std /** - * Computes the componentwise minimum of two - * vectorized data fields. The result is - * return as vectorized array in the form - * {min(x[0],y[0]), min(x[1],y[1]), - * ...}. + * Computes the componentwise minimum of two vectorized data fields. The + * result is return as vectorized array in the form {min(x[0],y[0]), + * min(x[1],y[1]), ...}. * * @relates VectorizedArray */ diff --git a/deal.II/include/deal.II/fe/mapping.h b/deal.II/include/deal.II/fe/mapping.h index 647de86ce8..4067d70983 100644 --- a/deal.II/include/deal.II/fe/mapping.h +++ b/deal.II/include/deal.II/fe/mapping.h @@ -591,7 +591,7 @@ public: mapping_contravariant_gradient() and mapping_piola_gradient() are only true as stated for linear mappings. If, for example, the mapping is bilinear then there is a missing - term associated with the derivative of of J. + term associated with the derivative of J. */ virtual void diff --git a/deal.II/include/deal.II/lac/parallel_block_vector.h b/deal.II/include/deal.II/lac/parallel_block_vector.h index 747d12a33e..3be362a635 100644 --- a/deal.II/include/deal.II/lac/parallel_block_vector.h +++ b/deal.II/include/deal.II/lac/parallel_block_vector.h @@ -226,6 +226,41 @@ namespace parallel void reinit (const BlockVector &V, const bool fast=false); + /** + * This function copies the data that has accumulated in the data buffer + * for ghost indices to the owning processor. For the meaning of the + * argument @p operation, see the entry on @ref GlossCompress + * "Compressing distributed vectors and matrices" in the glossary. + * + * There are two variants for this function. If called with argument @p + * VectorOperation::add adds all the data accumulated in ghost elements + * to the respective elements on the owning processor and clears the + * ghost array afterwards. If called with argument @p + * VectorOperation::insert, a set operation is performed. Since setting + * elements in a vector with ghost elements is ambiguous (as one can set + * both the element on the ghost site as well as the owning site), this + * operation makes the assumption that all data is set correctly on the + * owning processor. Upon call of compress(VectorOperation::insert), all + * ghost entries are therefore simply zeroed out (using + * zero_ghost_values()). In debug mode, a check is performed that makes + * sure that the data set is actually consistent between processors, + * i.e., whenever a non-zero ghost element is found, it is compared to + * the value on the owning processor and an exception is thrown if these + * elements do not agree. + * + */ + void compress (::dealii::VectorOperation::values operation); + + /** + * Fills the data field for ghost indices with the values stored in the + * respective positions of the owning processor. This function is needed + * before reading from ghosts. The function is @p const even though + * ghost data is changed. This is needed to allow functions with a @p + * const vector to perform the data exchange without creating + * temporaries. + */ + void update_ghost_values () const; + /** * Return whether the vector contains only elements with value * zero. This function is mainly for internal consistency checks and @@ -543,6 +578,29 @@ namespace parallel } + + template + inline + void + BlockVector::compress (::dealii::VectorOperation::values operation) + { + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).compress(operation); + } + + + + template + inline + void + BlockVector::update_ghost_values () const + { + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).update_ghost_values(); + } + + + template inline bool diff --git a/deal.II/include/deal.II/lac/sparsity_pattern.h b/deal.II/include/deal.II/lac/sparsity_pattern.h index 532162f06b..414e6ea92b 100644 --- a/deal.II/include/deal.II/lac/sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/sparsity_pattern.h @@ -842,7 +842,9 @@ public: /** * Copy data from an object of type CompressedSparsityPattern, - * CompressedSetSparsityPattern or CompressedSimpleSparsityPattern. + * CompressedSetSparsityPattern or CompressedSimpleSparsityPattern. Although + * not a compressed sparsity pattern, this function is also instantiated + * if the argument is of type SparsityPattern (i.e., the current class). * Previous content of this object is lost, and the sparsity pattern is in * compressed mode afterwards. */ diff --git a/deal.II/include/deal.II/numerics/data_out.h b/deal.II/include/deal.II/numerics/data_out.h index 7dfe9c9bd4..961fed7c20 100644 --- a/deal.II/include/deal.II/numerics/data_out.h +++ b/deal.II/include/deal.II/numerics/data_out.h @@ -253,7 +253,8 @@ private: void build_one_patch (const std::pair *cell_and_index, internal::DataOut::ParallelData &data, ::dealii::DataOutBase::Patch &patch, - const CurvedCellRegion curved_cell_region); + const CurvedCellRegion curved_cell_region, + std::vector > &patches); }; diff --git a/deal.II/include/deal.II/numerics/derivative_approximation.h b/deal.II/include/deal.II/numerics/derivative_approximation.h index 3fd49856c4..3b47dd3904 100644 --- a/deal.II/include/deal.II/numerics/derivative_approximation.h +++ b/deal.II/include/deal.II/numerics/derivative_approximation.h @@ -19,8 +19,11 @@ #include #include +#include +#include #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -673,12 +676,12 @@ private: template class DH, class InputVector, int spacedim> static void - approximate (const typename DH::active_cell_iterator &cell, + approximate (SynchronousIterators + ::active_cell_iterator,Vector::iterator> > const &cell, const Mapping &mapping, const DH &dof, const InputVector &solution, - const unsigned int component, - Vector &derivative_norm); + const unsigned int component); /** * Compute the derivative approximation on diff --git a/deal.II/source/base/data_out_base.cc b/deal.II/source/base/data_out_base.cc index d4d7a03852..e5518afcb8 100644 --- a/deal.II/source/base/data_out_base.cc +++ b/deal.II/source/base/data_out_base.cc @@ -6523,9 +6523,8 @@ void DataOutInterface::write_vtu_in_parallel (const char *filename MPI_File fh; err = MPI_File_open(comm, const_cast(filename), MPI_MODE_CREATE | MPI_MODE_WRONLY, info, &fh); - std::cout << err << std::endl; AssertThrow(err==0, ExcMessage("Unable to open file with MPI_File_open!")); - + MPI_File_set_size(fh, 0); // delete the file contents // this barrier is necessary, because otherwise others might already diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index df146a150b..669dccb71b 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -1115,16 +1115,17 @@ next_cell: + 1); - // next we have to build a mapping from the - // list of cells to their indices. for - // this, use the user_index field - std::vector saved_user_indices; - triangulation.save_user_indices (saved_user_indices); + // create a map pair -> SparsityPattern index + // TODO: we are no longer using user_indices for this because we can get + // pointer/index clashes when saving/restoring them. The following approach + // works, but this map can get quite big. Not sure about more efficient solutions. + std::map< std::pair, unsigned int > + indexmap; unsigned int index = 0; for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell, ++index) - cell->set_user_index (index); + indexmap[std::pair(cell->level(),cell->index())] = index; // next loop over all cells and // their neighbors to build the @@ -1153,21 +1154,15 @@ next_cell: && (cell->neighbor(f)->has_children() == false)) { - cell_connectivity.add (index, - cell->neighbor(f)->user_index()); - cell_connectivity.add (cell->neighbor(f)->user_index(), - index); + unsigned int other_index = indexmap.find( + std::pair(cell->neighbor(f)->level(),cell->neighbor(f)->index()))->second; + cell_connectivity.add (index, other_index); + cell_connectivity.add (other_index, index); } } - // now compress the so-built connectivity - // pattern and restore user indices. the - // const-cast is necessary since we treat - // the triangulation as constant (we here - // return it to its original state) + // now compress the so-built connectivity pattern cell_connectivity.compress (); - const_cast&>(triangulation) - .load_user_indices (saved_user_indices); } diff --git a/deal.II/source/numerics/data_out.cc b/deal.II/source/numerics/data_out.cc index 03e0cce91a..be256797bb 100644 --- a/deal.II/source/numerics/data_out.cc +++ b/deal.II/source/numerics/data_out.cc @@ -59,17 +59,11 @@ namespace internal /** - * In a WorkStream context, use this function to append the patch computed - * by the parallel stage to the array of patches. + * Dummy function used for WorkStream. */ template void - append_patch_to_list (const DataOutBase::Patch &patch, - std::vector > &patches) - { - patches.push_back (patch); - patches.back().patch_index = patches.size()-1; - } + copy(const DataOutBase::Patch &patch) {} } } @@ -81,7 +75,8 @@ DataOut:: build_one_patch (const std::pair *cell_and_index, internal::DataOut::ParallelData &data, DataOutBase::Patch &patch, - const CurvedCellRegion curved_cell_region) + const CurvedCellRegion curved_cell_region, + std::vector > &patches) { // use ucd_to_deal map as patch vertices are in the old, unnatural // ordering. if the mapping does not preserve locations @@ -306,6 +301,14 @@ build_one_patch (const std::pair *cell_and_index, patch.neighbors[f] = (*data.cell_to_patch_index_map)[neighbor->level()][neighbor->index()]; } + + const unsigned int patch_idx = cell_and_index->second; + // did we mess up the indices? + Assert(patch_idx < patches.size(), ExcInternalError()); + + // Put the patch in the patches vector + patches[patch_idx] = patch; + patches[patch_idx].patch_index = patch_idx; } @@ -357,19 +360,8 @@ void DataOut::build_patches (const Mapping > all_cells; { - // set the index of the first cell. if first_locally_owned_cell / - // next_locally_owned_cell returns non-active cells, then the index is not - // usable anyway, but otherwise we should keep track where we are - unsigned int index; cell_iterator cell = first_locally_owned_cell(); - if ((cell == this->triangulation->end()) - || - (cell->has_children())) - index = 0; - else - index = std::distance (this->triangulation->begin_active(), - active_cell_iterator(cell)); - for ( ; cell != this->triangulation->end(); ) + for (unsigned int index = 0; cell != this->triangulation->end(); ++index) { Assert (static_cast(cell->level()) < cell_to_patch_index_map.size(), @@ -381,25 +373,12 @@ void DataOut::build_patches (const Mappinglevel()][cell->index()] = all_cells.size(); all_cells.push_back (std::make_pair(cell, index)); - - // if both this and the next cell are active, then increment the index - // that keeps track on which active cell we are sitting correctly. if - // one of the cells is not active, then this index doesn't mean - // anything anyway, so just ignore it. same if we are at the end of - // the range - cell_iterator next = next_locally_owned_cell(cell); - if (!cell->has_children() && - next != this->triangulation->end() && - !next->has_children()) - index += std::distance (active_cell_iterator(cell), - active_cell_iterator(next)); - cell = next; + cell = next_locally_owned_cell(cell); } } this->patches.clear (); - this->patches.reserve (all_cells.size()); - Assert (this->patches.size() == 0, ExcInternalError()); + this->patches.resize(all_cells.size()); // now create a default object for the WorkStream object to work with unsigned int n_datasets=this->cell_data.size(); @@ -449,9 +428,9 @@ void DataOut::build_patches (const Mapping::build_one_patch, *this, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3, - curved_cell_region), - std_cxx1x::bind(&internal::DataOut::append_patch_to_list, - std_cxx1x::_1, std_cxx1x::ref(this->patches)), + curved_cell_region,std_cxx1x::ref(this->patches)), + std_cxx1x::bind(&internal::DataOut::copy, + std_cxx1x::_1), thread_data, sample_patch); } @@ -486,7 +465,7 @@ typename DataOut::cell_iterator DataOut::first_locally_owned_cell () { typename DataOut::cell_iterator - cell = this->triangulation->begin_active (); + cell = first_cell(); // skip cells if the current one has no children (is active) and is a ghost // or artificial cell diff --git a/deal.II/source/numerics/derivative_approximation.cc b/deal.II/source/numerics/derivative_approximation.cc index 397e688c60..cb8825c807 100644 --- a/deal.II/source/numerics/derivative_approximation.cc +++ b/deal.II/source/numerics/derivative_approximation.cc @@ -650,20 +650,25 @@ approximate_derivative (const Mapping &mapping, // Only act on the locally owned cells typedef FilteredIterator::active_cell_iterator> CellFilter; - - // There is no need for a copier because there is no conflict between threads + + typedef std_cxx1x::tuple::active_cell_iterator,Vector::iterator> + Iterators; + SynchronousIterators begin(Iterators (CellFilter(IteratorFilters::LocallyOwnedCell(), + dof_handler.begin_active()),derivative_norm.begin())), + end(Iterators (CellFilter(IteratorFilters::LocallyOwnedCell(),dof_handler.end()), + derivative_norm.end())); + + // There is no need for a copier because there is no conflict between threads // to write in derivative_norm. Scratch and CopyData are also useless. - WorkStream::run(CellFilter(IteratorFilters::LocallyOwnedCell(),dof_handler.begin_active()), - CellFilter(IteratorFilters::LocallyOwnedCell(),dof_handler.end()), - static_cast::active_cell_iterator const&, + WorkStream::run(begin,end, + static_cast const&, internal::Assembler::Scratch const&,internal::Assembler::CopyData &)> > (std_cxx1x::bind(DerivativeApproximation::template approximate, std_cxx1x::_1, std_cxx1x::cref(mapping), std_cxx1x::cref(dof_handler), - std_cxx1x::cref(solution),component, - std_cxx1x::ref(derivative_norm))), + std_cxx1x::cref(solution),component)), static_cast > (internal::Assembler::copier),internal::Assembler::Scratch (), internal::Assembler::CopyData ()); @@ -674,27 +679,21 @@ approximate_derivative (const Mapping &mapping, template class DH, class InputVector, int spacedim> void -DerivativeApproximation::approximate (const typename DH::active_cell_iterator &cell, +DerivativeApproximation::approximate (SynchronousIterators + ::active_cell_iterator,Vector::iterator> > const &cell, const Mapping &mapping, const DH &dof_handler, const InputVector &solution, - const unsigned int component, - Vector &derivative_norm) + const unsigned int component) { - // iterators over all cells and the - // respective entries in the output - // vector: - Vector::iterator derivative_norm_on_this_cell = derivative_norm.begin() + - std::distance(dof_handler.begin_active(),cell); - typename DerivativeDescription::Derivative derivative; // call the function doing the actual // work on this cell DerivativeApproximation::template approximate_cell - (mapping,dof_handler,solution,component,cell,derivative); + (mapping,dof_handler,solution,component,std_cxx1x::get<0>(cell.iterators),derivative); // evaluate the norm and fill the vector - *derivative_norm_on_this_cell - = DerivativeDescription::derivative_norm (derivative); + //*derivative_norm_on_this_cell + *std_cxx1x::get<1>(cell.iterators) = DerivativeDescription::derivative_norm (derivative); } diff --git a/deal.II/source/numerics/matrix_tools.cc b/deal.II/source/numerics/matrix_tools.cc index af0e940cc8..75d21fe9e8 100644 --- a/deal.II/source/numerics/matrix_tools.cc +++ b/deal.II/source/numerics/matrix_tools.cc @@ -16,7 +16,6 @@ #include #include -#include #include #include #include @@ -603,6 +602,41 @@ namespace MatrixCreator data.dof_indices, *matrix); } + + + + namespace AssemblerBoundary + { + struct Scratch + { + Scratch() {} + }; + + template + struct CopyData + { + CopyData() {}; + + CopyData(CopyData const &data); + + unsigned int dofs_per_cell; + std::vector dofs; + std::vector > dof_is_on_face; + typename DH::active_cell_iterator cell; + std::vector > cell_matrix; + std::vector > cell_vector; + }; + + template + CopyData::CopyData(CopyData const &data) : + dofs_per_cell(data.dofs_per_cell), + dofs(data.dofs), + dof_is_on_face(data.dof_is_on_face), + cell(data.cell), + cell_matrix(data.cell_matrix), + cell_vector(data.cell_vector) + {} + } } } @@ -834,37 +868,30 @@ namespace MatrixCreator { template void - create_boundary_mass_matrix_1 (std_cxx1x::tuple &, - const DoFHandler &, - const Quadrature & > commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const coefficient, - const std::vector &component_mapping, - const MatrixCreator::internal::IteratorRange > range, - Threads::Mutex &mutex) + create_boundary_mass_matrix_1 (typename DoFHandler::active_cell_iterator const &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch, + MatrixCreator::internal::AssemblerBoundary::CopyData > ©_data, + Mapping const &mapping, + FiniteElement const &fe, + Quadrature const &q, + typename FunctionMap::type const &boundary_functions, + Function const *const coefficient, + std::vector const &component_mapping) + { // All assertions for this function // are in the calling function // before creating threads. - const Mapping &mapping = std_cxx1x::get<0>(commons); - const DoFHandler &dof = std_cxx1x::get<1>(commons); - const Quadrature& q = std_cxx1x::get<2>(commons); - - const FiniteElement &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 cell_matrix(dofs_per_cell, dofs_per_cell); - Vector 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 | @@ -883,236 +910,246 @@ namespace MatrixCreator std::vector rhs_values_scalar (fe_values.n_quadrature_points); std::vector > rhs_values_system (fe_values.n_quadrature_points, - Vector(n_function_components)); - - std::vector dofs (dofs_per_cell); - std::vector dofs_on_face_vector (dofs_per_face); + Vector(n_function_components)); - // for each dof on the cell, have a - // flag whether it is on the face - std::vector dof_is_on_face(dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = range.first; - for (; cell!=range.second; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - // check if this face is on that part of - // the boundary we are interested in - if (boundary_functions.find(cell->face(face)->boundary_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 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::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 (copy_data.dofs_per_cell, + copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); + fe_values.reinit (cell, face); - // Special treatment - // for Hdiv and Hcurl - // elements, where only - // the normal or - // tangential component - // should be projected. - std::vector > normal_adjustment(fe_values.n_quadrature_points, - std::vector(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 &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 > normal_adjustment(fe_values.n_quadrature_points, + std::vector(n_components, 1.)); - if (!base.conforms(FiniteElementData::H1) && - base.conforms(FiniteElementData::Hdiv)) - for (unsigned int point=0; point &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; pointface(face)->boundary_indicator()) - ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); + if (!base.conforms(FiniteElementData::H1) && + base.conforms(FiniteElementData::Hdiv)) + for (unsigned int point=0; pointvalue_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointget_dof_indices (dofs); - cell->face(face)->get_dof_indices (dofs_on_face_vector); - for (unsigned int i=0; iface(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; pointface(face)->get_dof_indices (dofs_on_face_vector); + // for each dof on the cell, have a + // flag whether it is on the face + copy_data.dof_is_on_face.push_back(std::vector (copy_data.dofs_per_cell)); + // check for each of the dofs on this cell + // whether it is on the face + for (unsigned int i=0; i + void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData > const ©_data, + typename FunctionMap::type const &boundary_functions, + std::vector const &dof_to_boundary_mapping, + SparseMatrix &matrix, + Vector &rhs_vector) + { + // now transfer cell matrix and vector to the whole boundary matrix + // + // in the following: dof[i] holds the global index of the i-th degree of + // freedom on the present cell. If it is also a dof on the boundary, it + // must be a nonzero entry in the dof_to_boundary_mapping and then + // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. + // + // if dof[i] is not on the boundary, it should be zero on the boundary + // therefore on all quadrature points and finally all of its + // entries in the cell matrix and vector should be zero. If not, we + // throw an error (note: because of the evaluation of the shape + // functions only up to machine precision, the term "must be zero" + // really should mean: "should be very small". since this is only an + // assertion and not part of the code, we may choose "very small" + // quite arbitrarily) + // + // the main problem here is that the matrix or vector entry should also + // be zero if the degree of freedom dof[i] is on the boundary, but not + // on the present face, i.e. on another face of the same cell also + // on the boundary. We can therefore not rely on the + // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to + // determine whether dof[i] is a dof on the present face. We do so + // by getting the dofs on the face into @p{dofs_on_face_vector}, + // a vector as always. Usually, searching in a vector is + // inefficient, so we copy the dofs into a set, which enables binary + // searches. + unsigned int pos(0); + for (unsigned int face=0; face::faces_per_cell; ++face) + { + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) != + boundary_functions.end()) + { + for (unsigned int i=0; i void - create_boundary_mass_matrix_1<2,3> (std_cxx1x::tuple &, - const DoFHandler<2,3> &, - const Quadrature<1> & > , - SparseMatrix &, - const FunctionMap<3>::type &, - Vector &, - std::vector &, - const Function<3> *const , - const std::vector &, - const MatrixCreator::internal::IteratorRange > , - 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 > ©_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 const &component_mapping) { Assert(false,ExcNotImplemented()); } @@ -1121,17 +1158,17 @@ namespace MatrixCreator template <> void - create_boundary_mass_matrix_1<1,3> (std_cxx1x::tuple &, - const DoFHandler<1,3> &, - const Quadrature<0> & > , - SparseMatrix &, - const FunctionMap<3>::type &, - Vector &, - std::vector &, - const Function<3> *const , - const std::vector &, - const MatrixCreator::internal::IteratorRange > , - 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 > ©_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 const &component_mapping) { Assert(false,ExcNotImplemented()); } @@ -1184,47 +1221,24 @@ namespace MatrixCreator 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::active_cell_iterator active_cell_iterator; - std::vector > thread_ranges - = Threads::split_range (dof.begin_active(), - dof.end(), n_threads); - - // mutex to synchronise access to - // the matrix - Threads::Mutex mutex; - - typedef std_cxx1x::tuple &, - const DoFHandler &, - const Quadrature&> Commons; - - // then assemble in parallel - typedef void (*create_boundary_mass_matrix_1_t) - (Commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const coefficient, - const std::vector &component_mapping, - const MatrixCreator::internal::IteratorRange > range, - Threads::Mutex &mutex); - create_boundary_mass_matrix_1_t p - = &create_boundary_mass_matrix_1; - -//TODO: Use WorkStream here - for (unsigned int thread=0; thread > copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > + (std_cxx1x::bind(create_boundary_mass_matrix_1,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 > const &)> > (std_cxx1x::bind( + copy_boundary_mass_matrix_1,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); } @@ -1234,37 +1248,28 @@ namespace MatrixCreator template void - create_boundary_mass_matrix_1 (std_cxx1x::tuple &, - const hp::DoFHandler &, - const hp::QCollection &> commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const coefficient, - const std::vector &component_mapping, - const MatrixCreator::internal::IteratorRange > range, - Threads::Mutex &mutex) + create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler::active_cell_iterator const + &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch, + MatrixCreator::internal::AssemblerBoundary + ::CopyData > ©_data, + hp::MappingCollection const &mapping, + hp::FECollection const &fe_collection, + hp::QCollection const &q, + const typename FunctionMap::type &boundary_functions, + Function const *const coefficient, + std::vector const &component_mapping) { - const hp::MappingCollection &mapping = std_cxx1x::get<0>(commons); - const hp::DoFHandler &dof = std_cxx1x::get<1>(commons); - const hp::QCollection& q = std_cxx1x::get<2>(commons); - const hp::FECollection &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(0); - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != hp::DoFHandler::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - Assert (max_element == matrix.n()-1, ExcInternalError()); - }; -#endif + const FiniteElement &fe = cell->get_fe(); + const unsigned int dofs_per_face = fe.dofs_per_face; + + copy_data.cell = cell; + copy_data.dofs_per_cell = fe.dofs_per_cell; + copy_data.dofs.resize(copy_data.dofs_per_cell); + cell->get_dof_indices (copy_data.dofs); const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell(), max_dofs_per_face = fe_collection.max_dofs_per_face(); @@ -1287,148 +1292,131 @@ namespace MatrixCreator std::vector rhs_values_scalar; std::vector > rhs_values_system; - std::vector dofs (max_dofs_per_cell); - std::vector dofs_on_face_vector (max_dofs_per_face); + std::vector dofs_on_face_vector (dofs_per_face); - // for each dof on the cell, have a - // flag whether it is on the face - std::vector 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::active_cell_iterator cell = range.first; - for (; cell!=range.second; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - // check if this face is on that part of - // the boundary we are interested in - if (boundary_functions.find(cell->face(face)->boundary_indicator()) != - boundary_functions.end()) - { - x_fe_values.reinit (cell, face); - const FEFaceValues &fe_values = x_fe_values.get_present_fe_values (); + for (unsigned int face=0; face::faces_per_cell; ++face) + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(cell->face(face)->boundary_indicator()) != + boundary_functions.end()) + { + x_fe_values.reinit (cell, face); - const FiniteElement &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 &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 (copy_data.dofs_per_cell, + copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); - if (fe_is_system) - // FE has several components - { - rhs_values_system.resize (fe_values.n_quadrature_points, - Vector(n_function_components)); - boundary_functions.find(cell->face(face)->boundary_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(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(n_components)); - coefficient->vector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - for (unsigned int point=0; pointn_components==1) { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; ivalue_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; pointface(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(n_components)); + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_vector_values); + for (unsigned int point=0; pointvalue_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointface(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; pointface(face)->get_dof_indices (dofs_on_face_vector, + cell->active_fe_index()); + // for each dof on the cell, have a + // flag whether it is on the face + copy_data.dof_is_on_face.push_back(std::vector (copy_data.dofs_per_cell)); + // check for each of the dofs on this cell + // whether it is on the face + for (unsigned int i=0; iget_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 + void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary + ::CopyData > const ©_data, + typename FunctionMap::type const &boundary_functions, + std::vector const &dof_to_boundary_mapping, + SparseMatrix &matrix, + Vector &rhs_vector) + { + // now transfer cell matrix and vector to the whole boundary matrix + // + // in the following: dof[i] holds the global index of the i-th degree of + // freedom on the present cell. If it is also a dof on the boundary, it + // must be a nonzero entry in the dof_to_boundary_mapping and then + // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. + // + // if dof[i] is not on the boundary, it should be zero on the boundary + // therefore on all quadrature points and finally all of its + // entries in the cell matrix and vector should be zero. If not, we + // throw an error (note: because of the evaluation of the shape + // functions only up to machine precision, the term "must be zero" + // really should mean: "should be very small". since this is only an + // assertion and not part of the code, we may choose "very small" + // quite arbitrarily) + // + // the main problem here is that the matrix or vector entry should also + // be zero if the degree of freedom dof[i] is on the boundary, but not + // on the present face, i.e. on another face of the same cell also + // on the boundary. We can therefore not rely on the + // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to + // determine whether dof[i] is a dof on the present face. We do so + // by getting the dofs on the face into @p{dofs_on_face_vector}, + // a vector as always. Usually, searching in a vector is + // inefficient, so we copy the dofs into a set, which enables binary + // searches. + unsigned int pos(0); + for (unsigned int face=0; face::faces_per_cell; ++face) + { + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) != + boundary_functions.end()) + { #ifdef DEBUG - double max_diag_entry = 0; - for (unsigned int i=0; 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(0); + for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != hp::DoFHandler::invalid_dof_index) && + (*i > max_element)) + max_element = *i; + Assert (max_element == matrix.n()-1, ExcInternalError()); + + double max_diag_entry = 0; + for (unsigned int i=0; i max_diag_entry) + max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i)); #endif - // lock the matrix - Threads::Mutex::ScopedLock lock (mutex); - for (unsigned int i=0; i threads; - - // define starting and end point - // for each thread - typedef typename hp::DoFHandler::active_cell_iterator active_cell_iterator; - std::vector > thread_ranges - = Threads::split_range (dof.begin_active(), - dof.end(), n_threads); - - typedef std_cxx1x::tuple &, - const hp::DoFHandler &, - const hp::QCollection&> 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 &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const coefficient, - const std::vector &component_mapping, - const MatrixCreator::internal::IteratorRange > range, - Threads::Mutex &mutex); - create_boundary_mass_matrix_1_t p = &create_boundary_mass_matrix_1; - -//TODO: Use WorkStream here - for (unsigned int thread=0; thread > copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > + (std_cxx1x::bind(create_hp_boundary_mass_matrix_1,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 > const &)> > (std_cxx1x::bind( + copy_hp_boundary_mass_matrix_1,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); }