#include <grid/tria.h>
#include <grid/tria_iterator.h>
+#include <grid/intergrid_map.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_constraints.h>
}
+
+
+template <int dim>
+void
+DoFTools::compute_intergrid_constraints (const DoFHandler<dim> &coarse_grid,
+ const unsigned int coarse_component,
+ const DoFHandler<dim> &fine_grid,
+ const unsigned int fine_component,
+ const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+ ConstraintMatrix &constraints)
+{
+ // aliases to the finite elements
+ // used by the dof handlers:
+ const FiniteElement<dim> &coarse_fe = coarse_grid.get_fe(),
+ &fine_fe = fine_grid.get_fe();
+
+ // global numbers of dofs
+ const unsigned int n_coarse_dofs = coarse_grid.n_dofs(),
+ n_fine_dofs = fine_grid.n_dofs();
+
+ // local numbers of dofs
+ const unsigned int coarse_dofs_per_cell = coarse_fe.dofs_per_cell,
+ fine_dofs_per_cell = fine_fe.dofs_per_cell;
+
+ // alias the number of dofs per
+ // cell belonging to the
+ // coarse_component which is to be
+ // the restriction of the fine
+ // grid:
+ const unsigned int coarse_dofs_per_cell_component
+ = coarse_fe.base_element(coarse_fe.component_to_base(coarse_component)).dofs_per_cell;
+
+
+ // check whether component numbers
+ // are valid
+ Assert (coarse_component < coarse_fe.n_components(),
+ ExcInvalidComponent (coarse_component, coarse_fe.n_components()));
+ Assert (fine_component < fine_fe.n_components(),
+ ExcInvalidComponent (fine_component, fine_fe.n_components()));
+ // check whether respective finite
+ // elements are equal
+ Assert (coarse_fe.base_element (coarse_fe.component_to_base(coarse_component))
+ ==
+ fine_fe.base_element (fine_fe.component_to_base(fine_component)),
+ ExcFiniteElementsDontMatch());
+
+
+/*
+ * From here on: the term `parameter' refers to the selected component
+ * on the coarse grid and its analogon on the fine grid. The naming of
+ * variables containing this term is due to the fact that
+ * `selected_component' is longer, but also due to the fact that the
+ * code of this function was initially written for a program where the
+ * component which we wanted to match between grids was actually the
+ * `parameter' variable.
+ *
+ * Likewise, the terms `parameter grid' and `state grid' refer to the
+ * coarse and fine grids, respectively.
+ *
+ * Changing the names of variables would in principle be a good idea,
+ * but would not make things simpler and would be another source of
+ * errors. If anyone feels like doing so: patches would be welcome!
+ */
+
+
+
+ // set up vectors of cell-local
+ // data; each vector represents one
+ // degree of freedom of the
+ // coarse-grid variable in the
+ // fine-grid element
+ vector<Vector<double> > parameter_dofs (coarse_dofs_per_cell_component,
+ Vector<double>(fine_dofs_per_cell));
+ // for each coarse dof: find its
+ // position within the fine element
+ // and set this value to one in the
+ // respective vector (all other values
+ // are zero by construction)
+ for (unsigned int local_coarse_dof=0;
+ local_coarse_dof<coarse_dofs_per_cell_component;
+ ++local_coarse_dof)
+ parameter_dofs[local_coarse_dof]
+ (fine_fe.component_to_system_index (fine_component,local_coarse_dof)) = 1.;
+
+
+ // find out how many DoFs there are
+ // on the grids belonging to the
+ // components we want to match
+ unsigned int n_parameters_on_fine_grid=0;
+ if (true)
+ {
+ // have a flag for each dof on
+ // the fine grid and set it
+ // to true if this is an
+ // interesting dof. finally count
+ // how many true's there
+ vector<bool> dof_is_interesting (fine_grid.n_dofs(), false);
+ vector<unsigned int> local_dof_indices (fine_fe.dofs_per_cell);
+
+ for (DoFHandler<dim>::active_cell_iterator cell=fine_grid.begin_active();
+ cell!=fine_grid.end(); ++cell)
+ {
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
+ if (fine_fe.system_to_component_index(i).first == fine_component)
+ dof_is_interesting[local_dof_indices[i]] = true;
+ };
+
+ n_parameters_on_fine_grid = count (dof_is_interesting.begin(),
+ dof_is_interesting.end(),
+ true);
+ };
+
+ // vector to hold the representation of
+ // a single degree of freedom on the
+ // coarse grid (for the selected fe)
+ // on the fine grid
+ Vector<double> global_parameter_representation (n_fine_dofs);
+
+ // store the weights with which a dof
+ // on the parameter grid contributes
+ // to a dof on the fine grid. see the
+ // long doc below for more info
+ //
+ // allocate as many rows as there are
+ // parameter dofs on the coarse grid
+ // and as many columns as there are
+ // parameter dofs on the fine grid.
+ //
+ // weight_mapping is used to map the
+ // global (fine grid) parameter dof
+ // indices to the columns
+ FullMatrix<double> weights (n_coarse_dofs,
+ n_parameters_on_fine_grid);
+ // this is this mapping. there is one
+ // entry for each dof on the fine grid;
+ // if it is a parameter dof, then its
+ // value is the column in weights for
+ // that parameter dof, if it is any
+ // other dof, then its value is -1,
+ // indicating an error
+ vector<int> weight_mapping (n_fine_dofs, -1);
+
+ // set up this mapping
+ if (true)
+ {
+ // list of local dof numbers
+ // belonging to the parameter
+ // part of the FE on the fine grid
+ vector<unsigned int> parameter_dofs;
+ for (unsigned int i=0; i<coarse_dofs_per_cell; ++i)
+ parameter_dofs.push_back(fine_fe.component_to_system_index(coarse_component,i));
+
+ vector<unsigned int> local_dof_indices(fine_fe.dofs_per_cell);
+ unsigned int next_free_index=0;
+ for (typename DoFHandler<dim>::active_cell_iterator cell=fine_grid.begin_active();
+ cell != fine_grid.end(); ++cell)
+ {
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
+ // if this DoF is a
+ // parameter dof and has
+ // not yet been numbered,
+ // then do so
+ if ((fine_fe.system_to_component_index(i).first == fine_component) &&
+ (weight_mapping[local_dof_indices[i]] == -1))
+ {
+ weight_mapping[local_dof_indices[i]] = next_free_index;
+ ++next_free_index;
+ };
+ };
+
+ Assert (next_free_index == n_parameters_on_fine_grid,
+ ExcInternalError());
+ };
+
+
+ // for each cell on the parameter grid:
+ // find out which degrees of freedom on the
+ // fine grid correspond in which way to
+ // the degrees of freedom on the parameter
+ // grid
+ //
+ // since for continuous FEs some
+ // dofs exist on more than one
+ // cell, we have to track which
+ // ones were already visited. the
+ // problem is that if we visit a
+ // dof first on one cell and
+ // compute its weight with respect
+ // to some global dofs to be
+ // non-zero, and later visit the
+ // dof again on another cell and
+ // (since we are on another cell)
+ // recompute the weights with
+ // respect to the same dofs as
+ // above to be zero now, we have to
+ // preserve them. we therefore
+ // overwrite all weights if they
+ // are nonzero and do not enforce
+ // zero weights since that might be
+ // only due to the fact that we are
+ // on another cell.
+ //
+ // example:
+ // coarse grid
+ // | | |
+ // *-----*-----*
+ // | cell|cell |
+ // | 1 | 2 |
+ // | | |
+ // 0-----1-----*
+ //
+ // fine grid
+ // | | | | |
+ // *--*--*--*--*
+ // | | | | |
+ // *--*--*--*--*
+ // | | | | |
+ // *--x--y--*--*
+ //
+ // when on cell 1, we compute the
+ // weights of dof 'x' to be 1/2
+ // from parameter dofs 0 and 1,
+ // respectively. however, when
+ // later we are on cell 2, we again
+ // compute the prolongation of
+ // shape function 1 restricted to
+ // cell 2 to the globla grid and
+ // find that the weight of global
+ // dof 'x' now is zero. however, we
+ // should not overwrite the old
+ // value.
+ //
+ // we therefore always only set
+ // nonzero values. why adding up is
+ // not useful: dof 'y' would get
+ // weight 1 from parameter dof 1 on
+ // both cells 1 and 2, but the
+ // correct weight is nevertheless
+ // only 1.
+ typename DoFHandler<dim>::active_cell_iterator cell;
+ vector<unsigned int> parameter_dof_indices (coarse_fe.dofs_per_cell);
+
+ for (cell=coarse_grid.begin_active(); cell!=coarse_grid.end(); ++cell)
+ {
+ // get the global indices of the
+ // parameter dofs on this parameter
+ // grid cell
+ cell->get_dof_indices (parameter_dof_indices);
+
+ // loop over all dofs on this
+ // cell and check whether they
+ // are interesting for us
+ for (unsigned int local_parameter_dof=0;
+ local_parameter_dof<coarse_fe.dofs_per_cell;
+ ++local_parameter_dof)
+ if (coarse_fe.system_to_component_index(local_parameter_dof).first
+ ==
+ coarse_component)
+ {
+ global_parameter_representation.clear ();
+
+ // distribute the representation of
+ // #local_parameter_dof# on the
+ // parameter grid cell #cell# to
+ // the global data space
+ coarse_to_fine_grid_map[cell]->
+ set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof],
+ global_parameter_representation);
+ // now that we've got the global
+ // representation of each parameter
+ // dof, we've only got to clobber the
+ // non-zero entries in that vector and
+ // store the result
+ //
+ // what we have learned: if entry #i#
+ // of the global vector holds the value
+ // #v[i]#, then this is the weight with
+ // which the present dof contributes
+ // to #i#. there may be several such
+ // #i#s and their weights' sum should
+ // be one. Then, #v[i]# should be
+ // equal to #\sum_j w_{ij} p[j]# with
+ // #p[j]# be the values of the degrees
+ // of freedom on the coarse grid. we
+ // can thus compute constraints which
+ // link the degrees of freedom #v[i]#
+ // on the fine grid to those on the
+ // coarse grid, #p[j]#. Now to use
+ // these as real constraints, rather
+ // than as additional equations, we
+ // have to identify representants
+ // among the #i# for each #j#. this will
+ // be done by simply taking the first
+ // #i# for which #w_{ij}==1#.
+ for (unsigned int i=0; i<global_parameter_representation.size(); ++i)
+ // set this weight if it belongs
+ // to a parameter dof.
+ if (weight_mapping[i] != -1)
+ {
+ // only overwrite old
+ // value if not by
+ // zero
+ if (global_parameter_representation(i) != 0)
+ {
+ const unsigned int wi = parameter_dof_indices[local_parameter_dof],
+ wj = weight_mapping[i];
+ weights(wi,wj) = global_parameter_representation(i);
+ };
+ }
+ else
+ Assert (global_parameter_representation(i) == 0,
+ ExcInternalError());
+ };
+ };
+
+ // ok, now we have all weights for each
+ // dof on the fine grid. if in debug
+ // mode lets see if everything went smooth,
+ // i.e. each dof has sum of weights one
+ //
+ // in other words this means that
+ // if the sum of all shape
+ // functions on the parameter grid
+ // is one (which is always the
+ // case), then the representation
+ // on the state grid should be as
+ // well (division of unity)
+ //
+ // if the parameter grid has more
+ // than one component, then the
+ // respective dofs of the other
+ // components have sum of weights
+ // zero, of course. we do not
+ // explicitely ask which component
+ // a dof belongs to, but this at
+ // least tests some errors
+#ifdef DEBUG
+ for (unsigned int col=0; col<weights.n(); ++col)
+ {
+ double sum=0;
+ for (unsigned int row=0; row<weights.m(); ++row)
+ sum += weights(row,col);
+ Assert ((sum==1) ||
+ ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
+ };
+#endif
+
+
+ // now we know that the weights in each
+ // row constitute a constraint. enter this
+ // into the constraints object
+ //
+ // first task: for each parameter dof on
+ // the parameter grid, find a representant
+ // on the fine, global grid. this is
+ // possible since we use conforming
+ // finite element. we take this representant
+ // to be the first element in this row with
+ // weight identical to one
+ vector<int> representants(weights.m(), -1);
+ for (unsigned int parameter_dof=0; parameter_dof<weights.m(); ++parameter_dof)
+ {
+ unsigned int column=0;
+ for (; column<weights.n(); ++column)
+ if (weights(parameter_dof,column) == 1)
+ break;
+ Assert (column < weights.n(), ExcInternalError());
+
+ // now we know in which column of
+ // weights the representant is, but
+ // we don't know its global index. get
+ // it using the inverse operation of
+ // the weight_mapping
+ unsigned int global_dof=0;
+ for (; global_dof<weight_mapping.size(); ++global_dof)
+ if (weight_mapping[global_dof] == static_cast<int>(column))
+ break;
+ Assert (global_dof < weight_mapping.size(), ExcInternalError());
+
+ // now enter the representants global
+ // index into our list
+ representants[parameter_dof] = global_dof;
+ };
+
+ for (unsigned int global_dof=0; global_dof<n_fine_dofs; ++global_dof)
+ if (weight_mapping[global_dof] != -1)
+ // this global dof is a parameter
+ // dof, so it may carry a constraint
+ // note that for each global dof,
+ // the sum of weights shall be one,
+ // so we can find out whether this
+ // dof is constrained in the following
+ // way: if the only weight in this row
+ // is a one, and the representant for
+ // the parameter dof of the line in
+ // which this one is is the present
+ // dof, then we consider this dof
+ // to be unconstrained. otherwise,
+ // all other dofs are constrained
+ {
+ unsigned int first_used_row=0;
+ for (; first_used_row<weights.m(); ++first_used_row)
+ if (weights(first_used_row,weight_mapping[global_dof]) != 0)
+ break;
+
+ if ((weights(first_used_row,weight_mapping[global_dof]) == 1) &&
+ (representants[first_used_row] == static_cast<int>(global_dof)))
+ // dof unconstrained
+ continue;
+
+ // otherwise enter all constraints
+ constraints.add_line (global_dof);
+ for (unsigned int row=first_used_row; row<weights.m(); ++row)
+ if (weights(row,weight_mapping[global_dof]) != 0)
+ constraints.add_entry (global_dof,
+ representants[row],
+ weights(row,weight_mapping[global_dof]));
+ };
+};
+
+
+
+
+
+
// explicit instantiations
#if deal_II_dimension > 1
template void
const vector<bool>& local_select,
vector<bool>& selected_dofs);
+
+template
+void
+DoFTools::compute_intergrid_constraints (const DoFHandler<deal_II_dimension> &,
+ const unsigned int ,
+ const DoFHandler<deal_II_dimension> &,
+ const unsigned int ,
+ const InterGridMap<DoFHandler,deal_II_dimension> &,
+ ConstraintMatrix &);