* is the same as that for the alike triangulation iterators.
*
* This class also provides functions to create the sparsity patterns of
- * global matrices as well as matrices living on (parts of) the boundary
- * and handles some simple forms of transfer of data from one grid to another.
+ * global matrices as well as matrices living on (parts of) the boundary.
*
*
* \subsection{Distribution of indices for degrees of freedom}
* relevant in this context.
*
*
- * \subsection{Data transfer between grids}
- *
- * {\bf Note: The functionality described in this section has never been tested
- * and is expected to be broken or at least incomplete. Contact the author if
- * you need these functions for more information.}
- *
- * The #DoFHandler# class offers two functions #make_transfer_matrix# which create
- * a matrix to transform the data of one grid to another. The functions assumes the
- * coarsest mesh of the two grids to be the same. However there are few ways to
- * check this (only the number of cells on the coarsest grid is compared). Also,
- * the selected finite element type of the two degree of freedom handler objects
- * must be the same.
- *
- * The algorithm goes recursively from the coarse mesh cells to their children
- * until the grids differ at this level. It then tries to prolong or restrict the
- * old cell(s) to the new cell(s) and makes up a matrix of these prolongations and
- * restrictions. This matrix multiplied with a vector on the old grid yields an
- * approximation of the projection of the function on the old grid to the new one.
- *
- * Building and using the transfer matrix is usually quite an expensive operation,
- * since we have to perform two runs over all cells (one for building the sparsity
- * structure, one to build the entries) and because of the memory consumption.
- * It may, however, pay if you have many
- * equations, since then the entries in the matrix can be considered as block
- * entries which are then applied to all function values at a given degree of
- * freedom.
- *
- * To build the matrix, you first have to call
- * #make_transfer_matrix (old_dof_object, sparsity_pattern);#, then create a
- * sparse matrix out of this pattern, e.g. by #SparseMatrix<double> m(sparsity_pattern);#
- * and finally give this to the second run:
- * #make_transfer_matrix (old_dof_object, m);#. The spasity pattern created
- * by the first run is automatically compressed.
- *
- * When creating the #SparseMatrixStruct# sparsity pattern, you have to give the
- * dimension and the maximum number of entries per row. Obviously the image
- * dimension is the number of dofs on the new grid (you can get this using the
- * #n_dofs()# function), while the range dimension is the number of dofs on the
- * old grid. The maximum number of entries per row is determined by the maximum
- * number of levels $d$ which we have to cross upon transferring from one cell to
- * another (presently, transfer of one cell is only possible for #d=0,1#, i.e.
- * the two cells match or one is refined once more than the other, the
- * number of degrees of freedom per vertex $d_v$, those on lines $d_l$, those
- * on quads $d_q$ and the number of subcells a cell is
- * refined to, which is $2^{dim}$. The maximum number of entries per row in one
- * dimension is then given by $(2*d_l+d_v)*2+1$ if $d=1$. For example, a one
- * dimensional linear element would need two entries per row.
- * In two dimensions, the maxmimum number is $(4*d_q+12*d_l+5*d_v)*4+1$ if $d=1$.
- * You can get these numbers by drawing little pictures and counting, there is
- * no mystique behind this. You can also get the right number by calling the
- * #max_transfer_entries (max_level_difference)# function. The actual number
- * depends on the finite element selected and may be much less, especially in
- * higher dimensions.
- *
- *
- *
- * @author Wolfgang Bangerth, 1998 */
+ * @author Wolfgang Bangerth, 1998
+ */
template <int dim>
class DoFHandler : public Subscriptor,
public DoFDimensionInfo<dim>
SparseMatrixStruct &sparsity) const;
- /**
- * Make up the transfer matrix which
- * transforms the data vectors from one
- * triangulation to the present one.
- * You have to pass references to the old
- * dof handler object and to a matrix
- * sparsity object. This function therefore
- * only makes up the sparsity pattern.
- *
- * The given sparsity pattern is
- * compressed at the end.
- *
- * In the matrix, row indices belong to
- * new dof numbers, column indices to the
- * ones on the old grid. Therefore,
- * multiplying this matrix by a vector
- * of the old grid yields the vector on
- * the new one.
- *
- * For more details see the general
- * documentation for this class.
- */
- void make_transfer_matrix (const DoFHandler<dim> &transfer_from,
- SparseMatrixStruct &transfer_pattern) const;
-
- /**
- * Make up the transfer matrix which
- * transforms the data vectors from one
- * triangulation to the present one.
- * You have to pass references to the old
- * dof handler object and to a matrix
- * object. This function therefore
- * builds the matrix itself
- *
- * The matrix object should be
- * associated with the sparsity pattern
- * constructed by the other
- * #make_transfer_matrix# object.
- *
- * For more details see the general
- * documentation for this class.
- */
- void make_transfer_matrix (const DoFHandler<dim> &transfer_from,
- SparseMatrix<double> &transfer_matrix) const;
-
/**
* Return the maximum number of
* degrees of freedom a degree of freedom
*/
unsigned int max_couplings_between_boundary_dofs () const;
- /**
- * Return the maximum number of entries
- * a row in a transfer matrix may contain
- * if any two cells of which the dofs are
- * to be transferred differ in refinement
- * level at most by #max_level_diff#.
- * It is assumed that the finite element
- * selected by the last call to
- * #distribute_dofs# is used also for
- * the transfer process.
- */
- unsigned int max_transfer_entries (const unsigned int max_level_diff) const;
-
/**
* Take a vector of values which live on
* cells (e.g. an error per cell) and
/**
* Exception
*/
- DeclException0 (ExcOnlyOnelevelTransferImplemented);
- /**
- * Exception
- */
DeclException0 (ExcInvalidBoundaryIndicator);
/**
* Exception
unsigned int distribute_dofs_on_cell (active_cell_iterator &cell,
unsigned int next_free_dof);
- /**
- * Make up part of the sparsity pattern of
- * the transfer matrix by looking at the
- * two cells given.
- */
- void transfer_cell (const cell_iterator &old_cell,
- const cell_iterator &new_cell,
- SparseMatrixStruct &transfer_pattern) const;
-
- /**
- * Make up part of the transfer matrix by
- * looking at the two cells given.
- */
- void transfer_cell (const cell_iterator &old_cell,
- const cell_iterator &new_cell,
- SparseMatrix<double> &transfer_matrix) const;
-
/**
* Space to store the DoF numbers for the
* different levels. Analogous to the
-
-template <int dim>
-void DoFHandler<dim>::make_transfer_matrix (const DoFHandler<dim> &transfer_from,
- SparseMatrixStruct &transfer_pattern) const {
- Assert (tria->n_cells(0) == transfer_from.tria->n_cells(0),
- ExcGridsDoNotMatch());
- Assert (*selected_fe == *transfer_from.selected_fe,
- ExcGridsDoNotMatch());
- // assert for once at the beginning the
- // the matrices have the correct sizes
-#ifdef DEBUG
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- Assert (selected_fe->prolongate(c).m() == selected_fe->total_dofs,
- ExcMatrixHasWrongSize(selected_fe->prolongate(c).m()));
- Assert (selected_fe->prolongate(c).n() == selected_fe->total_dofs,
- ExcMatrixHasWrongSize(selected_fe->prolongate(c).n()));
- Assert (selected_fe->restrict(c).m() == selected_fe->total_dofs,
- ExcMatrixHasWrongSize(selected_fe->restrict(c).m()));
- Assert (selected_fe->restrict(c).n() == selected_fe->total_dofs,
- ExcMatrixHasWrongSize(selected_fe->restrict(c).n()));
- };
-#endif
-
- cell_iterator old_cell = transfer_from.begin(0),
- new_cell = begin(0);
- unsigned int n_coarse_cells = tria->n_cells(0);
-
- // first run: make up sparsity structure
- for (unsigned int j=0; j<n_coarse_cells; ++j, ++old_cell, ++new_cell)
- transfer_cell (old_cell, new_cell, transfer_pattern);
-
- transfer_pattern.compress();
-};
-
-
-
-template <int dim>
-void DoFHandler<dim>::make_transfer_matrix (const DoFHandler<dim> &transfer_from,
- SparseMatrix<double> &transfer_matrix) const {
- cell_iterator old_cell = transfer_from.begin(0),
- new_cell = begin(0);
- unsigned int n_coarse_cells = tria->n_cells(0);
-
- // second run: make up matrix entries
- for (unsigned int j=0; j<n_coarse_cells; ++j, ++old_cell, ++new_cell)
- transfer_cell (old_cell, new_cell, transfer_matrix);
-};
-
-
-
-template <int dim>
-void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterator &old_cell,
- const typename DoFHandler<dim>::cell_iterator &new_cell,
- SparseMatrixStruct &transfer_pattern) const {
- if (!new_cell->active() && !old_cell->active())
- // both cells have children; go deeper
- for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
- transfer_cell (old_cell->child(child), new_cell->child(child),
- transfer_pattern);
- else
- if (new_cell->active() && old_cell->active())
- // both cells have no children
- {
- vector<int> old_dofs(selected_fe->total_dofs);
- vector<int> new_dofs(selected_fe->total_dofs);
- old_cell->get_dof_indices (old_dofs);
- new_cell->get_dof_indices (new_dofs);
-
- // copy dofs one-by-one
- for (unsigned int j=0; j<old_dofs.size(); ++j)
- transfer_pattern.add (new_dofs[j], old_dofs[j]);
- }
- else
- if (!new_cell->active() && old_cell->active())
- // new cell has children, old one has not
- {
- cell_iterator child[GeometryInfo<dim>::children_per_cell];
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- child[c] = new_cell->child(c);
- Assert (child[c]->active(),
- ExcOnlyOnelevelTransferImplemented());
- };
-
- // numbers of old dofs
- vector<int> old_dof_indices (selected_fe->total_dofs);
- old_cell->get_dof_indices (old_dof_indices);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- // numbers of child dofs
- vector<int> child_dof_indices(selected_fe->total_dofs);
- child[c]->get_dof_indices (child_dof_indices);
-
- for (unsigned int k=0; k<selected_fe->total_dofs; ++k)
- for (unsigned int j=0; j<selected_fe->total_dofs; ++j)
- if (selected_fe->prolongate(c)(k,j) != 0.)
- transfer_pattern.add (child_dof_indices[k],
- old_dof_indices[j]);
- };
- } else {
- // old cell has children, new one has not
- cell_iterator child[GeometryInfo<dim>::children_per_cell];
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- child[c] = old_cell->child(c);
- Assert (child[c]->active(),
- ExcOnlyOnelevelTransferImplemented());
- };
-
- // numbers of new dofs
- vector<int> new_dof_indices(selected_fe->total_dofs);
- new_cell->get_dof_indices(new_dof_indices);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- // numbers of child dofs
- vector<int> child_dof_indices (selected_fe->total_dofs);
- child[c]->get_dof_indices (child_dof_indices);
-
- for (unsigned int k=0; k<selected_fe->total_dofs; ++k)
- for (unsigned int j=0; j<selected_fe->total_dofs; ++j)
- if (selected_fe->restrict(c)(k,j) != 0.)
- transfer_pattern.add (new_dof_indices[k],
- child_dof_indices[j]);
- };
- };
-};
-
-
-
-template <int dim>
-void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterator &old_cell,
- const typename DoFHandler<dim>::cell_iterator &new_cell,
- SparseMatrix<double> &transfer_matrix) const {
- if (!new_cell->active() && !old_cell->active())
- // both cells have children; go deeper
- for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
- transfer_cell (old_cell->child(child), new_cell->child(child),
- transfer_matrix);
- else
- if (new_cell->active() && old_cell->active())
- // both cells have no children
- {
- vector<int> old_dofs (selected_fe->total_dofs);
- vector<int> new_dofs (selected_fe->total_dofs);
- old_cell->get_dof_indices (old_dofs);
- new_cell->get_dof_indices (new_dofs);
-
- // copy dofs one-by-one
- for (unsigned int j=0; j<old_dofs.size(); ++j)
- // use the SparseMatrix<double>:: as a workaround
- // for a bug in egcs
- transfer_matrix.SparseMatrix<double>::set (new_dofs[j], old_dofs[j], 1.0);
- }
- else
- if (!new_cell->active() && old_cell->active())
- // new cell has children, old one has not
- {
- cell_iterator child[GeometryInfo<dim>::children_per_cell];
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- child[c] = new_cell->child(c);
- Assert (child[c]->active(),
- ExcOnlyOnelevelTransferImplemented());
- };
-
- // numbers of old dofs
- vector<int> old_dof_indices (selected_fe->total_dofs);
- old_cell->get_dof_indices (old_dof_indices);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- // numbers of child dofs
- vector<int> child_dof_indices (selected_fe->total_dofs);
- child[c]->get_dof_indices (child_dof_indices);
-
- for (unsigned int k=0; k<selected_fe->total_dofs; ++k)
- for (unsigned int j=0; j<selected_fe->total_dofs; ++j)
- if (selected_fe->prolongate(c)(k,j) != 0.)
- // use the SparseMatrix<double>::
- // as a workaround
- // for a bug in egcs
- transfer_matrix.SparseMatrix<double>::set (child_dof_indices[k],
- old_dof_indices[j],
- selected_fe->prolongate(c)(k,j));
- };
- } else {
- // old cell has children, new one has not
- cell_iterator child[GeometryInfo<dim>::children_per_cell];
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- child[c] = old_cell->child(c);
- Assert (child[c]->active(),
- ExcOnlyOnelevelTransferImplemented());
- };
-
- // numbers of new dofs
- vector<int> new_dof_indices (selected_fe->total_dofs);
- new_cell->get_dof_indices(new_dof_indices);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- {
- // numbers of child dofs
- vector<int> child_dof_indices (selected_fe->total_dofs);
- child[c]->get_dof_indices (child_dof_indices);
-
- for (unsigned int k=0; k<selected_fe->total_dofs; ++k)
- for (unsigned int j=0; j<selected_fe->total_dofs; ++j)
- if (selected_fe->restrict(c)(k,j) != 0.)
- // use the SparseMatrix<double>:: as
- // a workaround
- // for a bug in egcs
-
- transfer_matrix.SparseMatrix<double>::set (new_dof_indices[k],
- child_dof_indices[j],
- selected_fe->restrict(c)(k,j));
- };
- };
-};
-
-
-
#if deal_II_dimension == 1
template <>
return 0;
};
-
-
-template <>
-unsigned int DoFHandler<1>::max_transfer_entries (const unsigned int max_level_diff) const {
- Assert (max_level_diff<2, ExcOnlyOnelevelTransferImplemented());
- switch (max_level_diff)
- {
- case 0:
- return 1;
- case 1:
- return (2*selected_fe->dofs_per_line+
- selected_fe->dofs_per_vertex) * 2 + 1;
- };
- return 0;
-};
-
#endif
return 3*selected_fe->dofs_per_vertex + 2*selected_fe->dofs_per_line;
};
-
-
-template <>
-unsigned int DoFHandler<2>::max_transfer_entries (const unsigned int max_level_diff) const {
- Assert (max_level_diff<2, ExcOnlyOnelevelTransferImplemented());
- switch (max_level_diff)
- {
- case 0:
- return 1;
- case 1:
- return (4*selected_fe->dofs_per_quad +
- 12*selected_fe->dofs_per_line+
- 5*selected_fe->dofs_per_vertex) * 4 +1;
- };
- return 0;
-};
-
#endif