* 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
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<vector<bool> > &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<int> &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<int> &dof_to_boundary_mapping,
- SparseMatrixStruct &sparsity) const;
/**
/**
* Exception
*/
- DeclException2 (ExcDifferentDimensions,
- int, int,
- << "One dimension of the matrices is differing: "
- << arg1 << " vs. " << arg2);
- /**
- * Exception
- */
DeclException0 (ExcNoFESelected);
/**
* Exception
/**
* 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 << ") "
#include <vector>
-
-
+// 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.
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<int dim>
+ static void make_sparsity_pattern (const DoFHandler<dim>& 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<int dim>
+ static void make_sparsity_pattern (const DoFHandler<dim>& dof,
+ const vector<vector<bool> > &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<int dim>
+ static void
+ make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
+ const vector<int> &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<int dim>
+ static void
+ make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
+ const DoFHandler<dim>::FunctionMap &boundary_indicators,
+ const vector<int> &dof_to_boundary_mapping,
+ SparseMatrixStruct &sparsity);
+
/**
* Generate sparsity pattern for fluxes.
* This is a replacement of the
*/
template<int dim>
static void make_flux_sparsity_pattern (const DoFHandler<dim>&,
- SparseMatrixStruct &);
+ SparseMatrixStruct &);
/**
* Extract the indices of the degrees
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()));
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()));
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()));
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
#endif
-
-template <int dim>
-void DoFHandler<dim>::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<int> 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<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- }
-}
-
-template <int dim>
-void
-DoFHandler<dim>::make_sparsity_pattern (const vector<vector<bool> > &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; i<mask.size(); ++i)
- Assert (mask[i].size() == selected_fe->n_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<vector<bool> > dof_mask(dofs_per_cell,
- vector<bool>(dofs_per_cell, false));
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- dof_mask[i][j] = mask
- [selected_fe->system_to_component_index(i).first]
- [selected_fe->system_to_component_index(j).first];
-
-
- vector<int> 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<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (dof_mask[i][j] == true)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- };
-};
-
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-void DoFHandler<1>::make_boundary_sparsity_pattern (const vector<int> &,
- SparseMatrixStruct &) const {
- Assert (selected_fe != 0, ExcNoFESelected());
- Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void DoFHandler<1>::make_boundary_sparsity_pattern (const FunctionMap &,
- const vector<int> &,
- SparseMatrixStruct &) const {
- Assert (selected_fe != 0, ExcNoFESelected());
- Assert (false, ExcInternalError());
-};
-
-#endif
-
-
-
-template <int dim>
-void DoFHandler<dim>::make_boundary_sparsity_pattern (const vector<int> &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<int> 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<dofs_per_face; ++i)
- for (unsigned int j=0; j<dofs_per_face; ++j)
- sparsity.add (dof_to_boundary_mapping[dofs_on_this_face[i]],
- dof_to_boundary_mapping[dofs_on_this_face[j]]);
- };
-};
-
-
-
-template <int dim>
-void DoFHandler<dim>::make_boundary_sparsity_pattern (const FunctionMap &boundary_indicators,
- const vector<int> &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<int> 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<dofs_per_face; ++i)
- for (unsigned int j=0; j<dofs_per_face; ++j)
- sparsity.add (dof_to_boundary_mapping[dofs_on_this_face[i]],
- dof_to_boundary_mapping[dofs_on_this_face[j]]);
- };
-};
-
-
-
-
#if deal_II_dimension == 1
template <>
#include <grid/dof_constraints.h>
#include <fe/fe.h>
#include <numerics/dof_renumbering.h>
+#include <basic/dof_tools.h>
#include <lac/sparsematrix.h>
#include <vector>
// 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)
{
#include <basic/dof_tools.h>
#include <lac/sparsematrix.h>
+
+
+template <int dim>
+void
+DoFTools::make_sparsity_pattern (const DoFHandler<dim>& 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<int> dofs_on_this_cell(dofs_per_cell);
+ DoFHandler<dim>::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<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ }
+}
+
+
+template <int dim>
+void
+DoFTools::make_sparsity_pattern (const DoFHandler<dim>& dof,
+ const vector<vector<bool> > &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<mask.size(); ++i)
+ Assert (mask[i].size() == dof.get_fe().n_components(),
+ ExcDimensionMismatch(mask[i].size(), dof.get_fe().n_components()));
+
+ // first build a mask for each dof,
+ // not like the one given which represents
+ // components
+ vector<vector<bool> > dof_mask(dofs_per_cell,
+ vector<bool>(dofs_per_cell, false));
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ dof_mask[i][j] = mask
+ [dof.get_fe().system_to_component_index(i).first]
+ [dof.get_fe().system_to_component_index(j).first];
+
+
+ vector<int> dofs_on_this_cell(dofs_per_cell);
+ DoFHandler<dim>::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<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (dof_mask[i][j] == true)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ }
+}
+
+
#if deal_II_dimension == 1
-//TODO: Implement make_flux_sparsity_pattern in 1d
+template <>
+void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1>&,
+ const vector<int> &,
+ SparseMatrixStruct &)
+{
+ Assert (false, ExcInternalError());
+};
+
+
+template <>
+void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1>&,
+ const DoFHandler<1>::FunctionMap &,
+ const vector<int> &,
+ SparseMatrixStruct &)
+{
+ Assert (false, ExcInternalError());
+}
+
+#endif
+
+//TODO: Reinsert the last assertion: could not find max_element (G->W)
+
+template <int dim>
+void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
+ const vector<int> &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<int> dofs_on_this_face(dofs_per_face);
+ DoFHandler<dim>::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<dofs_per_face; ++i)
+ for (unsigned int j=0; j<dofs_per_face; ++j)
+ sparsity.add (dof_to_boundary_mapping[dofs_on_this_face[i]],
+ dof_to_boundary_mapping[dofs_on_this_face[j]]);
+ };
+};
+
+
+
+template <int dim>
+void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
+ const DoFHandler<dim>::FunctionMap &boundary_indicators,
+ const vector<int> &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<dim>::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<int> dofs_on_this_face(dofs_per_face);
+ DoFHandler<dim>::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<dofs_per_face; ++i)
+ for (unsigned int j=0; j<dofs_per_face; ++j)
+ sparsity.add (dof_to_boundary_mapping[dofs_on_this_face[i]],
+ dof_to_boundary_mapping[dofs_on_this_face[j]]);
+ };
+};
+
+
+
+
+#if deal_II_dimension == 1
+
+//TODO: Implement make_flux_sparsity_pattern in 1D ?
template <>
void
DoFTools::make_flux_sparsity_pattern (const DoFHandler<1>&,
// explicit instantiations
#if deal_II_dimension > 1
-template void DoFTools::make_flux_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
- SparseMatrixStruct &sparsity);
+template void
+DoFTools::make_flux_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+ SparseMatrixStruct &sparsity);
+template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+ const vector<int> &,
+ SparseMatrixStruct &);
+template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+ const DoFHandler<deal_II_dimension>::FunctionMap &boundary_indicators,
+ const vector<int> &dof_to_boundary_mapping,
+ SparseMatrixStruct &sparsity);
#endif
+
+template void
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+ SparseMatrixStruct &sparsity);
+
+template void
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+ const vector<vector<bool> > &mask,
+ SparseMatrixStruct &sparsity);
+
template void DoFTools::extract_dofs(const DoFHandler<deal_II_dimension>& dof,
const vector<bool>& local_select,
vector<bool>& selected_dofs);
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<int> dofs_on_this_cell(dofs_per_cell);
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<int> dofs_on_this_cell(dofs_per_cell);
#include <grid/dof_constraints.h>
#include <grid/tria_iterator.h>
#include <basic/data_out.h>
+#include <basic/dof_tools.h>
#include <base/function.h>
#include <fe/fe.h>
#include <base/quadrature.h>
template <int dim>
void ProblemBase<dim>::set_tria_and_dof (Triangulation<dim> *t,
- DoFHandler<dim> *d) {
+ DoFHandler<dim> *d)
+{
tria = t;
dof_handler = d;
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
#include <numerics/assembler.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
+#include <basic/dof_tools.h>
#include <lac/vector.h>
#include <lac/sparsematrix.h>
#include <lac/precondition.h>
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<double> mass_matrix (sparsity);
// 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