From 070b450338f1bb0f58033b5cb367e013b6afd59a Mon Sep 17 00:00:00 2001 From: kanschat Date: Fri, 3 Dec 1999 23:29:38 +0000 Subject: [PATCH] make_*_sparsity_pattern moved to DoFTools git-svn-id: https://svn.dealii.org/trunk@1975 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/dofs/dof_handler.h | 151 ------------ deal.II/deal.II/include/dofs/dof_tools.h | 181 ++++++++++++++- deal.II/deal.II/source/dofs/dof_handler.cc | 188 +-------------- .../deal.II/source/dofs/dof_renumbering.cc | 3 +- deal.II/deal.II/source/dofs/dof_tools.cc | 215 +++++++++++++++++- deal.II/deal.II/source/dofs/mg_dof_handler.cc | 4 +- .../source/multigrid/mg_dof_handler.cc | 4 +- deal.II/deal.II/source/numerics/base.cc | 6 +- deal.II/deal.II/source/numerics/vectors.cc | 9 +- 9 files changed, 412 insertions(+), 349 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 6bf15cb273..ae357b1e3d 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -202,28 +202,6 @@ class DoFDimensionInfo<3> { * into new degree of freedom indices. * * - * \subsection{Setting up sparsity patterns} - * - * When assembling system matrices, the entries are usually of the form - * $a_{ij} = a(\phi_i, \phi_j)$, where $a$ is a bilinear functional, often an - * integral. When using sparse matrices, we therefore only need to reserve space - * for those $a_{ij}$ only, which are nonzero, which is the same as to say that - * that the basis functions $\phi_i$ and $\phi_j$ have a nonempty intersection of - * their support. Since the support of basis functions is bound only on cells - * on which they are located or to which they are adjacent, to determine the - * sparsity pattern it is sufficient to loop over all cells and connect all - * basis functions on each cell with all other basis functions on that cell. - * There may be finite elements for which not all basis functions on a cell - * connect with each other, but no use of this case is made since no examples - * where this occurs are known to the author. - * - * When setting up sparsity patterns for matrices on the boundary, the same - * procedure is done, except for the fact that the loop only goes over faces - * on the boundary and the basis functions thereon. It is assumed that all - * other basis functions on a cell adjacent to the boundary vanish at the - * boundary itself, except for those which are located on the boundary. - * - * * \subsection{Boundaries} * * When projecting the traces of functions to the boundary or parts thereof, one @@ -432,121 +410,6 @@ class DoFHandler : public Subscriptor, void make_hanging_node_constraints (ConstraintMatrix &) const; - /** - * Write the sparsity structure of the - * full matrix including constrained - * degrees of freedom into the - * matrix structure. The sparsity pattern - * does not include entries introduced by - * the elimination of constrained nodes. - * The sparsity - * pattern is not compressed, since if - * you want to call - * #ConstraintMatrix::condense(1)# - * afterwards, new entries have to be - * added. However, if you don't want to call - * #ConstraintMatrix::condense(1)#, you - * have to compress the matrix yourself, - * using #SparseMatrixStruct::compress()#. - */ - void make_sparsity_pattern (SparseMatrixStruct &) const; - - /** - * This function does mistly the same as - * the above one, but it is specialized for - * mixed finite elements and allows to - * specify which variables couple in which - * equation. For example, if wanted to solve - * the Stokes equations, - * \begin{verbatim} - * -\Delta \vec u + \nabla p = 0, - * \div u = 0 - * \end{verbatim} - * in, say, two space dimensions, using - * nonstabilized Q2/Q2/Q1 mixed elements - * (i.e. using the #FESystem# class), then - * you don't want all degrees of freedom - * to couple in each equation. You rather - * may want to give the following pattern - * of couplings: - * \begin{verbatim} - * 1 0 1 - * 0 1 1 - * 1 1 0 - * \end{verbatim} - * where "1" indicates that two variables - * (i.e. components of the #FESystem#) - * couple in the respective equation, and - * a "0" means no coupling, in which case - * it is not necessary to allocate space - * in the matrix structure. Obviously, the - * mask refers to components of the - * composed #FESystem#, rather than to the - * degrees of freedom contained in there. - * - * This function is designed to accept - * a mask, like the one shown above, - * through the #mask# parameter, which - * contains boolean values. It builds - * the matrix structure just like the - * previous function, but does not create - * elements if not specified by the mask. - */ - void make_sparsity_pattern (const vector > &mask, - SparseMatrixStruct &) const; - - /** - * Write the sparsity structure of the - * matrix composed of the basis functions - * on the boundary into the - * matrix structure. The sparsity pattern - * does not include entries introduced by - * the elimination of constrained nodes. - * The sparsity - * pattern is not compressed, since if - * you want to call - * #ConstraintMatrix::condense(1)# - * afterwards, new entries have to be - * added. However, if you don't want to call - * #ConstraintMatrix::condense(1)#, you - * have to compress the matrix yourself, - * using #SparseMatrixStruct::compress()#. - * - * Since this function is obviously useless - * in one spatial dimension, it is not - * implemented. - */ - void make_boundary_sparsity_pattern (const vector &dof_to_boundary_mapping, - SparseMatrixStruct &) const; - - /** - * Write the sparsity structure of the - * matrix composed of the basis functions - * on the boundary into the - * matrix structure. In contrast to the - * previous function, only those parts - * of the boundary are considered of which - * the boundary indicator is listed in the - * set of numbers passed to this function. - * - * In fact, rather than a #set# of boundary - * indicators, a #map# needs to be passed, - * since most of the functions handling with - * boundary indicators take a mapping of - * boundary indicators and the respective - * boundary functions. The boundary function, - * however, is ignored in this function. - * If you have no functions at hand, but only - * the boundary indicators, set the function - * pointers to null pointers. - * - * Since this function is obviously useless - * in one spatial dimension, it is not - * implemented. - */ - void make_boundary_sparsity_pattern (const FunctionMap &boundary_indicators, - const vector &dof_to_boundary_mapping, - SparseMatrixStruct &sparsity) const; /** @@ -1270,13 +1133,6 @@ class DoFHandler : public Subscriptor, /** * Exception */ - DeclException2 (ExcDifferentDimensions, - int, int, - << "One dimension of the matrices is differing: " - << arg1 << " vs. " << arg2); - /** - * Exception - */ DeclException0 (ExcNoFESelected); /** * Exception @@ -1321,13 +1177,6 @@ class DoFHandler : public Subscriptor, /** * Exception */ - DeclException2 (ExcInvalidMaskDimension, - int, int, - << "The dimension of the mask " << arg1 << " does not match" - << " the number of components in the finite element object."); - /** - * Exception - */ DeclException2 (ExcInvalidComponent, int, int, << "The component you gave (" << arg1 << ") " diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index 5a0e3afb56..98eb87ca8a 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -8,8 +8,33 @@ #include - - +// This is part of the original documentation of DoFHandler +// I tried to put the essential part into the documentation of +// the function, where it belongs. Still, I keep it for a while, +// to be able to recover missing information. +/* + * \subsection{Setting up sparsity patterns} + * + * When assembling system matrices, the entries are usually of the form + * $a_{ij} = a(\phi_i, \phi_j)$, where $a$ is a bilinear functional, often an + * integral. When using sparse matrices, we therefore only need to reserve space + * for those $a_{ij}$ only, which are nonzero, which is the same as to say that + * that the basis functions $\phi_i$ and $\phi_j$ have a nonempty intersection of + * their support. Since the support of basis functions is bound only on cells + * on which they are located or to which they are adjacent, to determine the + * sparsity pattern it is sufficient to loop over all cells and connect all + * basis functions on each cell with all other basis functions on that cell. + * There may be finite elements for which not all basis functions on a cell + * connect with each other, but no use of this case is made since no examples + * where this occurs are known to the author. + * + * When setting up sparsity patterns for matrices on the boundary, the same + * procedure is done, except for the fact that the loop only goes over faces + * on the boundary and the basis functions thereon. It is assumed that all + * other basis functions on a cell adjacent to the boundary vanish at the + * boundary itself, except for those which are located on the boundary. + * + */ /** * Operations on DoF-numbers. @@ -24,6 +49,156 @@ class DoFTools { public: + /** + * Locate non-zero entries of the system matrix. + * + * This function computes the possible + * positions of non-zero entries + * in the global system + * matrix. We assume that a + * finite element basis function + * is non-zero on a cell only if + * its dual degree of freedom is + * associated with the interior, + * a face, an edge or a vertex of + * this cell. As a result, the + * matrix entry between two basis + * functions can be non-zero only + * if they correspond to degrees + * of freedom of at least one + * common + * cell. Therefore, + * #make_sparsity_pattern# just + * loops over all cells and + * enters all couplings local to + * that cell. + * + * Since this process is purely + * local, the sparsity pattern + * does not provide for entries introduced by + * the elimination of hanging nodes. + * They have to be taken care of + * by a call to + * #ConstraintMatrix::condense()# + * afterwards. + * + * Remember using + * #SparseMatrixStruct::compress()# + * after generating the pattern. + */ + template + static void make_sparsity_pattern (const DoFHandler& dof, + SparseMatrixStruct &); + + /** + * Locate non-zero entries for + * mixed methods. + * This function does mostly the same as + * the other + * #make_sparsity_pattern#, but + * it is specialized for + * mixed finite elements and allows to + * specify which variables couple in which + * equation. For example, if wanted to solve + * the Stokes equations, + * \begin{verbatim} + * -\Delta \vec u + \nabla p = 0, + * \div u = 0 + * \end{verbatim} + * in, two space dimensions, using + * stable Q2/Q1 mixed elements + * (using the #FESystem# class), then + * you don't want all degrees of freedom + * to couple in each equation. You rather + * may want to give the following pattern + * of couplings: + * \begin{verbatim} + * 1 0 1 + * 0 1 1 + * 1 1 0 + * \end{verbatim} + * where "1" indicates that two variables + * (i.e. components of the #FESystem#) + * couple in the respective equation, and + * a "0" means no coupling, in which case + * it is not necessary to allocate space + * in the matrix structure. Obviously, the + * mask refers to components of the + * composed #FESystem#, rather than to the + * degrees of freedom contained in there. + * + * This function is designed to accept + * a mask, like the one shown above, + * through the #mask# parameter, which + * contains boolean values. It builds + * the matrix structure just like the + * previous function, but does not create + * elements if not specified by the mask. + */ + template + static void make_sparsity_pattern (const DoFHandler& dof, + const vector > &mask, + SparseMatrixStruct &); + + /** + * Write the sparsity structure of the + * matrix composed of the basis functions + * on the boundary into the + * matrix structure. The sparsity pattern + * does not include entries introduced by + * the elimination of constrained nodes. + * The sparsity + * pattern is not compressed, since if + * you want to call + * #ConstraintMatrix::condense(1)# + * afterwards, new entries have to be + * added. However, if you don't want to call + * #ConstraintMatrix::condense(1)#, you + * have to compress the matrix yourself, + * using #SparseMatrixStruct::compress()#. + * + * Since this function is obviously useless + * in one spatial dimension, it is not + * implemented. + */ + template + static void + make_boundary_sparsity_pattern (const DoFHandler& dof, + const vector &dof_to_boundary_mapping, + SparseMatrixStruct &); + + /** + * Write the sparsity structure of the + * matrix composed of the basis functions + * on the boundary into the + * matrix structure. In contrast to the + * previous function, only those parts + * of the boundary are considered of which + * the boundary indicator is listed in the + * set of numbers passed to this function. + * + * In fact, rather than a #set# of boundary + * indicators, a #map# needs to be passed, + * since most of the functions handling with + * boundary indicators take a mapping of + * boundary indicators and the respective + * boundary functions. The boundary function, + * however, is ignored in this function. + * If you have no functions at hand, but only + * the boundary indicators, set the function + * pointers to null pointers. + * + * Since this function is obviously useless + * in one spatial dimension, it is not + * implemented. + */ + template + static void + make_boundary_sparsity_pattern (const DoFHandler& dof, + const DoFHandler::FunctionMap &boundary_indicators, + const vector &dof_to_boundary_mapping, + SparseMatrixStruct &sparsity); + /** * Generate sparsity pattern for fluxes. * This is a replacement of the @@ -38,7 +213,7 @@ class DoFTools */ template static void make_flux_sparsity_pattern (const DoFHandler&, - SparseMatrixStruct &); + SparseMatrixStruct &); /** * Extract the indices of the degrees diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index b7b9589a09..5501a2749f 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -1779,12 +1779,12 @@ void DoFHandler<2>::make_hanging_node_constraints (ConstraintMatrix &constraints Assert(2*selected_fe->dofs_per_vertex+selected_fe->dofs_per_line == selected_fe->constraints().n(), - ExcDifferentDimensions(2*selected_fe->dofs_per_vertex+ + ExcDimensionMismatch(2*selected_fe->dofs_per_vertex+ selected_fe->dofs_per_line, selected_fe->constraints().n())); Assert(selected_fe->dofs_per_vertex+2*selected_fe->dofs_per_line == selected_fe->constraints().m(), - ExcDifferentDimensions(3*selected_fe->dofs_per_vertex+ + ExcDimensionMismatch(3*selected_fe->dofs_per_vertex+ 2*selected_fe->dofs_per_line, selected_fe->constraints().m())); @@ -1872,7 +1872,7 @@ void DoFHandler<3>::make_hanging_node_constraints (ConstraintMatrix &constraints selected_fe->dofs_per_quad == selected_fe->constraints().n(), - ExcDifferentDimensions(4*selected_fe->dofs_per_vertex+ + ExcDimensionMismatch(4*selected_fe->dofs_per_vertex+ 4*selected_fe->dofs_per_line+ selected_fe->dofs_per_quad, selected_fe->constraints().n())); @@ -1881,7 +1881,7 @@ void DoFHandler<3>::make_hanging_node_constraints (ConstraintMatrix &constraints 4*selected_fe->dofs_per_quad == selected_fe->constraints().m(), - ExcDifferentDimensions(5*selected_fe->dofs_per_vertex+ + ExcDimensionMismatch(5*selected_fe->dofs_per_vertex+ 12*selected_fe->dofs_per_line+ 4*selected_fe->dofs_per_quad, selected_fe->constraints().m())); @@ -1949,11 +1949,11 @@ void DoFHandler<3>::make_hanging_node_constraints (ConstraintMatrix &constraints Assert (dofs_on_children.size() == selected_fe->constraints().m(), - ExcDifferentDimensions(dofs_on_children.size(), + ExcDimensionMismatch(dofs_on_children.size(), selected_fe->constraints().m())); Assert (dofs_on_mother.size() == selected_fe->constraints().n(), - ExcDifferentDimensions(dofs_on_mother.size(), + ExcDimensionMismatch(dofs_on_mother.size(), selected_fe->constraints().n())); // for each row in the constraint @@ -1972,182 +1972,6 @@ void DoFHandler<3>::make_hanging_node_constraints (ConstraintMatrix &constraints #endif - -template -void DoFHandler::make_sparsity_pattern (SparseMatrixStruct &sparsity) const -{ - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (sparsity.n_rows() == n_dofs(), - ExcDifferentDimensions (sparsity.n_rows(), n_dofs())); - Assert (sparsity.n_cols() == n_dofs(), - ExcDifferentDimensions (sparsity.n_cols(), n_dofs())); - - const unsigned int dofs_per_cell = selected_fe->dofs_per_cell; - vector dofs_on_this_cell(dofs_per_cell); - active_cell_iterator cell = begin_active(), - endc = end(); - for (; cell!=endc; ++cell) - { - cell->get_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this cell - for (unsigned int i=0; i -void -DoFHandler::make_sparsity_pattern (const vector > &mask, - SparseMatrixStruct &sparsity) const -{ - const unsigned int dofs_per_cell = selected_fe->dofs_per_cell; - - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (sparsity.n_rows() == n_dofs(), - ExcDifferentDimensions (sparsity.n_rows(), n_dofs())); - Assert (sparsity.n_cols() == n_dofs(), - ExcDifferentDimensions (sparsity.n_cols(), n_dofs())); - Assert (mask.size() == selected_fe->n_components(), - ExcInvalidMaskDimension(mask.size(), selected_fe->n_components())); - for (unsigned int i=0; in_components(), - ExcInvalidMaskDimension(mask[i].size(), selected_fe->n_components())); - - // first build a mask for each dof, - // not like the one given which represents - // components - vector > dof_mask(dofs_per_cell, - vector(dofs_per_cell, false)); - for (unsigned int i=0; isystem_to_component_index(i).first] - [selected_fe->system_to_component_index(j).first]; - - - vector dofs_on_this_cell(dofs_per_cell); - active_cell_iterator cell = begin_active(), - endc = end(); - for (; cell!=endc; ++cell) - { - cell->get_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this cell - for (unsigned int i=0; i -void DoFHandler<1>::make_boundary_sparsity_pattern (const vector &, - SparseMatrixStruct &) const { - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (false, ExcInternalError()); -}; - - - -template <> -void DoFHandler<1>::make_boundary_sparsity_pattern (const FunctionMap &, - const vector &, - SparseMatrixStruct &) const { - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (false, ExcInternalError()); -}; - -#endif - - - -template -void DoFHandler::make_boundary_sparsity_pattern (const vector &dof_to_boundary_mapping, - SparseMatrixStruct &sparsity) const { - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (dof_to_boundary_mapping.size() == n_dofs(), ExcInternalError()); - Assert (sparsity.n_rows() == n_boundary_dofs(), - ExcDifferentDimensions (sparsity.n_rows(), n_boundary_dofs())); - Assert (sparsity.n_cols() == n_boundary_dofs(), - ExcDifferentDimensions (sparsity.n_cols(), n_boundary_dofs())); - Assert (*max_element(dof_to_boundary_mapping.begin(), - dof_to_boundary_mapping.end()) == (signed int)sparsity.n_rows()-1, - ExcInternalError()); - - const unsigned int dofs_per_face = selected_fe->dofs_per_face; - vector dofs_on_this_face(dofs_per_face); - active_face_iterator face = begin_active_face(), - endf = end_face(); - for (; face!=endf; ++face) - if (face->at_boundary()) - { - face->get_dof_indices (dofs_on_this_face); - - // make sure all dof indices have a - // boundary index - Assert (*min_element(dofs_on_this_face.begin(), - dofs_on_this_face.end()) >=0, - ExcInternalError()); - - // make sparsity pattern for this cell - for (unsigned int i=0; i -void DoFHandler::make_boundary_sparsity_pattern (const FunctionMap &boundary_indicators, - const vector &dof_to_boundary_mapping, - SparseMatrixStruct &sparsity) const { - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (dof_to_boundary_mapping.size() == n_dofs(), ExcInternalError()); - Assert (boundary_indicators.find(255) == boundary_indicators.end(), - ExcInvalidBoundaryIndicator()); - Assert (sparsity.n_rows() == n_boundary_dofs(boundary_indicators), - ExcDifferentDimensions (sparsity.n_rows(), n_boundary_dofs(boundary_indicators))); - Assert (sparsity.n_cols() == n_boundary_dofs(boundary_indicators), - ExcDifferentDimensions (sparsity.n_cols(), n_boundary_dofs(boundary_indicators))); - Assert (*max_element(dof_to_boundary_mapping.begin(), - dof_to_boundary_mapping.end()) == (signed int)sparsity.n_rows()-1, - ExcInternalError()); - - const unsigned int dofs_per_face = selected_fe->dofs_per_face; - vector dofs_on_this_face(dofs_per_face); - active_face_iterator face = begin_active_face(), - endf = end_face(); - for (; face!=endf; ++face) - if (boundary_indicators.find(face->boundary_indicator()) != - boundary_indicators.end()) - { - face->get_dof_indices (dofs_on_this_face); - - // make sure all dof indices have a - // boundary index - Assert (*min_element(dofs_on_this_face.begin(), - dofs_on_this_face.end()) >=0, - ExcInternalError()); - // make sparsity pattern for this cell - for (unsigned int i=0; i diff --git a/deal.II/deal.II/source/dofs/dof_renumbering.cc b/deal.II/deal.II/source/dofs/dof_renumbering.cc index d2606a949b..6261c537bc 100644 --- a/deal.II/deal.II/source/dofs/dof_renumbering.cc +++ b/deal.II/deal.II/source/dofs/dof_renumbering.cc @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -23,7 +24,7 @@ void DoFRenumbering::Cuthill_McKee (DoFHandler &dof_handler, // make the connection graph SparseMatrixStruct sparsity (dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); - dof_handler.make_sparsity_pattern (sparsity); + DoFTools::make_sparsity_pattern (dof_handler, sparsity); if (use_constraints) { diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 0fa8bed12f..a6f4d88662 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -11,10 +11,199 @@ #include #include + + +template +void +DoFTools::make_sparsity_pattern (const DoFHandler& dof, + SparseMatrixStruct &sparsity) +{ + const unsigned int n_dofs = dof.n_dofs(); + + Assert (sparsity.n_rows() == n_dofs, + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + Assert (sparsity.n_cols() == n_dofs, + ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); + + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + vector dofs_on_this_cell(dofs_per_cell); + DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + cell->get_dof_indices (dofs_on_this_cell); + // make sparsity pattern for this cell + for (unsigned int i=0; i +void +DoFTools::make_sparsity_pattern (const DoFHandler& dof, + const vector > &mask, + SparseMatrixStruct &sparsity) +{ + const unsigned int n_dofs = dof.n_dofs(); + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + + Assert (sparsity.n_rows() == n_dofs, + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + Assert (sparsity.n_cols() == n_dofs, + ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); + Assert (mask.size() == dof.get_fe().n_components(), + ExcDimensionMismatch(mask.size(), dof.get_fe().n_components())); + for (unsigned int i=0; i > dof_mask(dofs_per_cell, + vector(dofs_per_cell, false)); + for (unsigned int i=0; i dofs_on_this_cell(dofs_per_cell); + DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + cell->get_dof_indices (dofs_on_this_cell); + // make sparsity pattern for this cell + for (unsigned int i=0; i +void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1>&, + const vector &, + SparseMatrixStruct &) +{ + Assert (false, ExcInternalError()); +}; + + +template <> +void +DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1>&, + const DoFHandler<1>::FunctionMap &, + const vector &, + SparseMatrixStruct &) +{ + Assert (false, ExcInternalError()); +} + +#endif + +//TODO: Reinsert the last assertion: could not find max_element (G->W) + +template +void +DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, + const vector &dof_to_boundary_mapping, + SparseMatrixStruct &sparsity) +{ + const unsigned int n_dofs = dof.n_dofs(); + + Assert (dof_to_boundary_mapping.size() == n_dofs, ExcInternalError()); + Assert (sparsity.n_rows() == dof.n_boundary_dofs(), + ExcDimensionMismatch (sparsity.n_rows(), dof.n_boundary_dofs())); + Assert (sparsity.n_cols() == dof.n_boundary_dofs(), + ExcDimensionMismatch (sparsity.n_cols(), dof.n_boundary_dofs())); +// Assert (*max_element(dof_to_boundary_mapping.begin(), +// dof_to_boundary_mapping.end()) == (signed int)sparsity.n_rows()-1, +// ExcInternalError()); + + const unsigned int dofs_per_face = dof.get_fe().dofs_per_face; + vector dofs_on_this_face(dofs_per_face); + DoFHandler::active_face_iterator face = dof.begin_active_face(), + endf = dof.end_face(); + for (; face!=endf; ++face) + if (face->at_boundary()) + { + face->get_dof_indices (dofs_on_this_face); + + // make sure all dof indices have a + // boundary index +// Assert (*min_element(dofs_on_this_face.begin(), +// dofs_on_this_face.end()) >=0, +// ExcInternalError()); + + // make sparsity pattern for this cell + for (unsigned int i=0; i +void DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, + const DoFHandler::FunctionMap &boundary_indicators, + const vector &dof_to_boundary_mapping, + SparseMatrixStruct &sparsity) +{ + const unsigned int n_dofs = dof.n_dofs(); + + Assert (dof_to_boundary_mapping.size() == n_dofs, ExcInternalError()); + Assert (boundary_indicators.find(255) == boundary_indicators.end(), + DoFHandler::ExcInvalidBoundaryIndicator()); + Assert (sparsity.n_rows() == dof.n_boundary_dofs(boundary_indicators), + ExcDimensionMismatch (sparsity.n_rows(), dof.n_boundary_dofs(boundary_indicators))); + Assert (sparsity.n_cols() == dof.n_boundary_dofs(boundary_indicators), + ExcDimensionMismatch (sparsity.n_cols(), dof.n_boundary_dofs(boundary_indicators))); +// Assert (*max_element(dof_to_boundary_mapping.begin(), +// dof_to_boundary_mapping.end()) == (signed int)sparsity.n_rows()-1, +// ExcInternalError()); + + const unsigned int dofs_per_face = dof.get_fe().dofs_per_face; + vector dofs_on_this_face(dofs_per_face); + DoFHandler::active_face_iterator face = dof.begin_active_face(), + endf = dof.end_face(); + for (; face!=endf; ++face) + if (boundary_indicators.find(face->boundary_indicator()) != + boundary_indicators.end()) + { + face->get_dof_indices (dofs_on_this_face); + + // make sure all dof indices have a + // boundary index +// Assert (*min_element(dofs_on_this_face.begin(), +// dofs_on_this_face.end()) >=0, +// ExcInternalError()); + // make sparsity pattern for this cell + for (unsigned int i=0; i void DoFTools::make_flux_sparsity_pattern (const DoFHandler<1>&, @@ -148,9 +337,29 @@ DoFTools::extract_level_dofs(const unsigned int level, // explicit instantiations #if deal_II_dimension > 1 -template void DoFTools::make_flux_sparsity_pattern (const DoFHandler& dof, - SparseMatrixStruct &sparsity); +template void +DoFTools::make_flux_sparsity_pattern (const DoFHandler& dof, + SparseMatrixStruct &sparsity); +template void +DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, + const vector &, + SparseMatrixStruct &); +template void +DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, + const DoFHandler::FunctionMap &boundary_indicators, + const vector &dof_to_boundary_mapping, + SparseMatrixStruct &sparsity); #endif + +template void +DoFTools::make_sparsity_pattern (const DoFHandler& dof, + SparseMatrixStruct &sparsity); + +template void +DoFTools::make_sparsity_pattern (const DoFHandler& dof, + const vector > &mask, + SparseMatrixStruct &sparsity); + template void DoFTools::extract_dofs(const DoFHandler& dof, const vector& local_select, vector& selected_dofs); diff --git a/deal.II/deal.II/source/dofs/mg_dof_handler.cc b/deal.II/deal.II/source/dofs/mg_dof_handler.cc index 0b4e5c4b48..499d92e29e 100644 --- a/deal.II/deal.II/source/dofs/mg_dof_handler.cc +++ b/deal.II/deal.II/source/dofs/mg_dof_handler.cc @@ -1544,9 +1544,9 @@ void MGDoFHandler::make_sparsity_pattern (const unsigned int level, SparseMatrixStruct &sparsity) const { Assert (selected_fe != 0, ExcNoFESelected()); Assert (sparsity.n_rows() == n_dofs(level), - ExcDifferentDimensions (sparsity.n_rows(), n_dofs(level))); + ExcDimensionMismatch (sparsity.n_rows(), n_dofs(level))); Assert (sparsity.n_cols() == n_dofs(level), - ExcDifferentDimensions (sparsity.n_cols(), n_dofs(level))); + ExcDimensionMismatch (sparsity.n_cols(), n_dofs(level))); const unsigned int dofs_per_cell = selected_fe->dofs_per_cell; vector dofs_on_this_cell(dofs_per_cell); diff --git a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc index 0b4e5c4b48..499d92e29e 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc @@ -1544,9 +1544,9 @@ void MGDoFHandler::make_sparsity_pattern (const unsigned int level, SparseMatrixStruct &sparsity) const { Assert (selected_fe != 0, ExcNoFESelected()); Assert (sparsity.n_rows() == n_dofs(level), - ExcDifferentDimensions (sparsity.n_rows(), n_dofs(level))); + ExcDimensionMismatch (sparsity.n_rows(), n_dofs(level))); Assert (sparsity.n_cols() == n_dofs(level), - ExcDifferentDimensions (sparsity.n_cols(), n_dofs(level))); + ExcDimensionMismatch (sparsity.n_cols(), n_dofs(level))); const unsigned int dofs_per_cell = selected_fe->dofs_per_cell; vector dofs_on_this_cell(dofs_per_cell); diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 53f7fe8fa2..8a0e99c71d 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -36,7 +37,8 @@ ProblemBase::ProblemBase () : template void ProblemBase::set_tria_and_dof (Triangulation *t, - DoFHandler *d) { + DoFHandler *d) +{ tria = t; dof_handler = d; @@ -88,7 +90,7 @@ void ProblemBase::assemble (const Equation &equation, constraints.clear (); dof_handler->make_hanging_node_constraints (constraints); constraints.close (); - dof_handler->make_sparsity_pattern (system_sparsity); + DoFTools::make_sparsity_pattern (*dof_handler, system_sparsity); constraints.condense (system_sparsity); // reinite system matrix diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index b28c76cf17..621168aadb 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -386,7 +387,7 @@ void VectorTools::project (const DoFHandler &dof, SparseMatrixStruct sparsity(dof.n_dofs(), dof.n_dofs(), dof.max_couplings_between_dofs()); - dof.make_sparsity_pattern (sparsity); + DoFTools::make_sparsity_pattern (dof, sparsity); constraints.condense (sparsity); SparseMatrix mass_matrix (sparsity); @@ -645,8 +646,10 @@ VectorTools::project_boundary_values (const DoFHandler &dof, // set up sparsity structure SparseMatrixStruct sparsity(dof.n_boundary_dofs(boundary_functions), dof.max_couplings_between_boundary_dofs()); - dof.make_boundary_sparsity_pattern (boundary_functions, dof_to_boundary_mapping, - sparsity); + DoFTools::make_boundary_sparsity_pattern (dof, + boundary_functions, + dof_to_boundary_mapping, + sparsity); // note: for three or more dimensions, there // may be constrained nodes on the boundary -- 2.39.5