const DoFHandler<dim> &fine_grid,
const unsigned int fine_component,
const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
- ConstraintMatrix &constraints,
- std::vector<std::map<unsigned int, float> > *transfer_representation)
+ ConstraintMatrix &constraints)
+{
+ // 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
+ //
+ // in the original implementation,
+ // the weights array was actually
+ // of FullMatrix<double> type. this
+ // wasted huge amounts of memory,
+ // but was fast. nonetheless, since
+ // the memory consumption was
+ // quadratic in the number of
+ // degrees of freedom, this was not
+ // very practical, so we now use a
+ // vector of rows of the matrix,
+ // and in each row a vector of
+ // pairs (colnum,value). this seems
+ // like the best tradeoff between
+ // memory and speed, as it is now
+ // linear in memory and still fast
+ // enough.
+ //
+ // to save some memory and since
+ // the weights are usually
+ // (negative) powers of 2, we
+ // choose the value type of the
+ // matrix to be @p{float} rather
+ // than @p{double}.
+ std::vector<std::map<unsigned int, float> > weights;
+
+ // 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
+ std::vector<int> weight_mapping;
+
+ const unsigned int n_parameters_on_fine_grid
+ = compute_intergrid_weights_1 (coarse_grid, coarse_component, fine_grid, fine_component,
+ coarse_to_fine_grid_map, weights, weight_mapping);
+
+ // global numbers of dofs
+ const unsigned int n_coarse_dofs = coarse_grid.n_dofs(),
+ n_fine_dofs = fine_grid.n_dofs();
+
+
+ // get an array in which we store
+ // which dof on the coarse grid is
+ // a parameter and which is not
+ std::vector<bool> coarse_dof_is_parameter (coarse_grid.n_dofs());
+ if (true)
+ {
+ std::vector<bool> mask (coarse_grid.get_fe().n_components(),
+ false);
+ mask[coarse_component] = true;
+ extract_dofs (coarse_grid, mask, coarse_dof_is_parameter);
+ };
+
+ // 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. the
+ // representant will become an
+ // unconstrained degree of freedom,
+ // while all others will be
+ // constrained to this dof (and
+ // possibly others)
+ std::vector<int> representants(n_coarse_dofs, -1);
+ for (unsigned int parameter_dof=0; parameter_dof<n_coarse_dofs;
+ ++parameter_dof)
+ if (coarse_dof_is_parameter[parameter_dof] == true)
+ {
+ // if this is the line of a
+ // parameter dof on the
+ // coarse grid, then it
+ // should have at least one
+ // dependent node on the fine
+ // grid
+ Assert (weights[parameter_dof].size() > 0, ExcInternalError());
+
+ // find the column where the
+ // representant is mentioned
+ std::map<unsigned int,float>::const_iterator i = weights[parameter_dof].begin();
+ for (; i!=weights[parameter_dof].end(); ++i)
+ if (i->second == 1)
+ break;
+ Assert (i!=weights[parameter_dof].end(), ExcInternalError());
+ const unsigned int column = i->first;
+
+ // 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;
+ }
+ else
+ {
+ // consistency check: if this
+ // is no parameter dof on the
+ // coarse grid, then the
+ // respective row must be
+ // empty!
+ Assert (weights[parameter_dof].size() == 0, ExcInternalError());
+ };
+
+
+
+ // note for people that want to
+ // optimize this function: the
+ // largest part of the computing
+ // time is spent in the following,
+ // rather innocent block of
+ // code. basically, it must be the
+ // ConstraintMatrix::add_entry call
+ // which takes the bulk of the
+ // time, but it is not known to the
+ // author how to make it faster...
+ std::vector<std::pair<unsigned int,double> > constraint_line;
+ 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
+ {
+ const unsigned int col = weight_mapping[global_dof];
+ Assert (col < n_parameters_on_fine_grid, ExcInternalError());
+
+ unsigned int first_used_row=0;
+ if (true)
+ {
+ std::map<unsigned int,float>::const_iterator col_entry;
+ for (; first_used_row<n_coarse_dofs; ++first_used_row)
+ {
+ col_entry = weights[first_used_row].find(col);
+ if (col_entry != weights[first_used_row].end())
+ break;
+ };
+
+ if ((col_entry->second == 1) &&
+ (representants[first_used_row] == static_cast<int>(global_dof)))
+ // dof unconstrained or
+ // constrained to itself
+ // (in case this cell is
+ // mapped to itself, rather
+ // than to children of
+ // itself)
+ continue;
+ };
+
+
+ // otherwise enter all constraints
+ constraints.add_line (global_dof);
+
+ constraint_line.clear ();
+ for (unsigned int row=first_used_row; row<n_coarse_dofs; ++row)
+ {
+ const std::map<unsigned int,float>::const_iterator
+ j = weights[row].find(col);
+ if ((j != weights[row].end()) && (j->second != 0))
+ constraint_line.push_back (std::make_pair(representants[row],
+ j->second));
+ };
+
+ constraints.add_entries (global_dof, constraint_line);
+ };
+};
+
+
+
+template <int dim>
+void
+DoFTools::
+compute_intergrid_transfer_representation (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,
+ std::vector<std::map<unsigned int, float> > &transfer_representation)
+{
+
+
+ // 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
+ //
+ // in the original implementation,
+ // the weights array was actually
+ // of FullMatrix<double> type. this
+ // wasted huge amounts of memory,
+ // but was fast. nonetheless, since
+ // the memory consumption was
+ // quadratic in the number of
+ // degrees of freedom, this was not
+ // very practical, so we now use a
+ // vector of rows of the matrix,
+ // and in each row a vector of
+ // pairs (colnum,value). this seems
+ // like the best tradeoff between
+ // memory and speed, as it is now
+ // linear in memory and still fast
+ // enough.
+ //
+ // to save some memory and since
+ // the weights are usually
+ // (negative) powers of 2, we
+ // choose the value type of the
+ // matrix to be @p{float} rather
+ // than @p{double}.
+ std::vector<std::map<unsigned int, float> > weights;
+
+ // 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
+ std::vector<int> weight_mapping;
+
+ compute_intergrid_weights_1 (coarse_grid, coarse_component, fine_grid, fine_component,
+ coarse_to_fine_grid_map, weights, weight_mapping);
+
+ // now compute the requested
+ // representation
+ transfer_representation.clear ();
+ transfer_representation.resize (weights.size());
+
+ const unsigned int n_global_parm_dofs
+ = std::count_if (weight_mapping.begin(), weight_mapping.end(),
+ std::bind2nd (std::not_equal_to<int> (), -1));
+
+ // first construct the inverse
+ // mapping of weight_mapping
+ std::vector<unsigned int> inverse_weight_mapping (n_global_parm_dofs,
+ DoFHandler<dim>::invalid_dof_index);
+ for (unsigned int i=0; i<weight_mapping.size(); ++i)
+ {
+ const unsigned int parameter_dof = weight_mapping[i];
+ // if this global dof is a
+ // parameter
+ if (parameter_dof != static_cast<unsigned int>(-1))
+ {
+ Assert (parameter_dof < n_global_parm_dofs, ExcInternalError());
+ Assert (inverse_weight_mapping[parameter_dof] == DoFHandler<dim>::invalid_dof_index,
+ ExcInternalError());
+
+ inverse_weight_mapping[parameter_dof] = i;
+ };
+ };
+
+ // next copy over weights array
+ // and replace respective
+ // numbers
+ const unsigned int n_coarse_dofs = coarse_grid.n_dofs();
+ for (unsigned int i=0; i<n_coarse_dofs; ++i)
+ {
+ std::map<unsigned int, float>::const_iterator j = weights[i].begin();
+ for (; j!=weights[i].end(); ++j)
+ transfer_representation[i][inverse_weight_mapping[j->first]] = j->second;
+ };
+};
+
+
+
+template <int dim>
+unsigned int
+DoFTools::compute_intergrid_weights_1 (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,
+ std::vector<std::map<unsigned int, float> > &weights,
+ std::vector<int> &weight_mapping)
{
// aliases to the finite elements
// used by the dof handlers:
dof_is_interesting[local_dof_indices[i]] = true;
};
- n_parameters_on_fine_grid = count (dof_is_interesting.begin(),
- dof_is_interesting.end(),
- true);
- };
-
- // get an array in which we store
- // which dof on the coarse grid is
- // a parameter and which is not
- std::vector<bool> coarse_dof_is_parameter (coarse_grid.n_dofs());
- if (true)
- {
- std::vector<bool> mask (coarse_grid.get_fe().n_components(),
- false);
- mask[coarse_component] = true;
- extract_dofs (coarse_grid, mask, coarse_dof_is_parameter);
- };
+ n_parameters_on_fine_grid = std::count (dof_is_interesting.begin(),
+ dof_is_interesting.end(),
+ true);
+ };
-
-
- // 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
- //
- // in the original implementation,
- // the weights array was actually
- // of FullMatrix<double> type. this
- // wasted huge amounts of memory,
- // but was fast. nonetheless, since
- // the memory consumption was
- // quadratic in the number of
- // degrees of freedom, this was not
- // very practical, so we now use a
- // vector of rows of the matrix,
- // and in each row a vector of
- // pairs (colnum,value). this seems
- // like the best tradeoff between
- // memory and speed, as it is now
- // linear in memory and still fast
- // enough.
- //
- // to save some memory and since
- // the weights are usually
- // (negative) powers of 2, we
- // choose the value type of the
- // matrix to be @p{float} rather
- // than @p{double}.
- std::vector<std::map<unsigned int, float> > weights(n_coarse_dofs);
- // 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
- std::vector<int> weight_mapping (n_fine_dofs, -1);
+ // set up the weights mapping
+ weights.clear ();
+ weights.resize (n_coarse_dofs);
- // set up this mapping
+ weight_mapping.clear ();
+ weight_mapping.resize (n_fine_dofs, -1);
+
if (true)
{
std::vector<unsigned int> local_dof_indices(fine_fe.dofs_per_cell);
// you want to read more
// information on the algorithm
// used.
- compute_intergrid_weights (coarse_grid, coarse_component,
- coarse_to_fine_grid_map, parameter_dofs,
- weight_mapping, weights);
+ compute_intergrid_weights_2 (coarse_grid, coarse_component,
+ coarse_to_fine_grid_map, parameter_dofs,
+ weight_mapping, weights);
+
// ok, now we have all weights for each
// dof on the fine grid. if in debug
};
#endif
- // if the user wants to have a
- // representation of the transfer
- // matrix, the provide it
- if (transfer_representation != 0)
- {
- transfer_representation->clear ();
- transfer_representation->resize (weights.size());
-
- const unsigned int n_global_parm_dofs
- = std::count_if (weight_mapping.begin(), weight_mapping.end(),
- std::bind2nd (std::not_equal_to<int> (), -1));
-
- // first construct the inverse
- // mapping of weight_mapping
- std::vector<unsigned int> inverse_weight_mapping (n_global_parm_dofs,
- DoFHandler<dim>::invalid_dof_index);
- for (unsigned int i=0; i<weight_mapping.size(); ++i)
- {
- const unsigned int parameter_dof = weight_mapping[i];
- // if this global dof is a
- // parameter
- if (parameter_dof != static_cast<unsigned int>(-1))
- {
- Assert (parameter_dof < n_global_parm_dofs, ExcInternalError());
- Assert (inverse_weight_mapping[parameter_dof] == DoFHandler<dim>::invalid_dof_index,
- ExcInternalError());
-
- inverse_weight_mapping[parameter_dof] = i;
- };
- };
-
- // next copy over weights array
- // and replace respective
- // numbers
- for (unsigned int i=0; i<n_coarse_dofs; ++i)
- {
- std::map<unsigned int, float>::const_iterator j = weights[i].begin();
- for (; j!=weights[i].end(); ++j)
- (*transfer_representation)[i][inverse_weight_mapping[j->first]] = j->second;
- };
- };
-
- // 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. the
- // representant will become an
- // unconstrained degree of freedom,
- // while all others will be
- // constrained to this dof (and
- // possibly others)
- std::vector<int> representants(n_coarse_dofs, -1);
- for (unsigned int parameter_dof=0; parameter_dof<n_coarse_dofs;
- ++parameter_dof)
- if (coarse_dof_is_parameter[parameter_dof] == true)
- {
- // if this is the line of a
- // parameter dof on the
- // coarse grid, then it
- // should have at least one
- // dependent node on the fine
- // grid
- Assert (weights[parameter_dof].size() > 0, ExcInternalError());
-
- // find the column where the
- // representant is mentioned
- std::map<unsigned int,float>::const_iterator i = weights[parameter_dof].begin();
- for (; i!=weights[parameter_dof].end(); ++i)
- if (i->second == 1)
- break;
- Assert (i!=weights[parameter_dof].end(), ExcInternalError());
- const unsigned int column = i->first;
-
- // 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;
- }
- else
- {
- // consistency check: if this
- // is no parameter dof on the
- // coarse grid, then the
- // respective row must be
- // empty!
- Assert (weights[parameter_dof].size() == 0, ExcInternalError());
- };
-
-
-
- // note for people that want to
- // optimize this function: the
- // largest part of the computing
- // time is spent in the following,
- // rather innocent block of
- // code. basically, it must be the
- // ConstraintMatrix::add_entry call
- // which takes the bulk of the
- // time, but it is not known to the
- // author how to make it faster...
- std::vector<std::pair<unsigned int,double> > constraint_line;
- 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
- {
- const unsigned int col = weight_mapping[global_dof];
- Assert (col < n_parameters_on_fine_grid, ExcInternalError());
-
- unsigned int first_used_row=0;
- if (true)
- {
- std::map<unsigned int,float>::const_iterator col_entry;
- for (; first_used_row<n_coarse_dofs; ++first_used_row)
- {
- col_entry = weights[first_used_row].find(col);
- if (col_entry != weights[first_used_row].end())
- break;
- };
-
- if ((col_entry->second == 1) &&
- (representants[first_used_row] == static_cast<int>(global_dof)))
- // dof unconstrained or
- // constrained to itself
- // (in case this cell is
- // mapped to itself, rather
- // than to children of
- // itself)
- continue;
- };
-
-
- // otherwise enter all constraints
- constraints.add_line (global_dof);
-
- constraint_line.clear ();
- for (unsigned int row=first_used_row; row<n_coarse_dofs; ++row)
- {
- const std::map<unsigned int,float>::const_iterator
- j = weights[row].find(col);
- if ((j != weights[row].end()) && (j->second != 0))
- constraint_line.push_back (std::make_pair(representants[row],
- j->second));
- };
-
- constraints.add_entries (global_dof, constraint_line);
- };
+ return n_parameters_on_fine_grid;
};
+
template <int dim>
void
-DoFTools::compute_intergrid_weights (const DoFHandler<dim> &coarse_grid,
- const unsigned int coarse_component,
- const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
- const std::vector<Vector<double> > ¶meter_dofs,
- const std::vector<int> &weight_mapping,
- std::vector<std::map<unsigned int,float> > &weights)
+DoFTools::compute_intergrid_weights_2 (const DoFHandler<dim> &coarse_grid,
+ const unsigned int coarse_component,
+ const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+ const std::vector<Vector<double> > ¶meter_dofs,
+ const std::vector<int> &weight_mapping,
+ std::vector<std::map<unsigned int,float> > &weights)
{
// simply distribute the range of
// cells to different threads
Threads::ThreadManager thread_manager;
for (unsigned int i=0; i<multithread_info.n_default_threads; ++i)
Threads::spawn (thread_manager,
- Threads::encapsulate (&DoFTools::template compute_intergrid_weights_1<dim>)
+ Threads::encapsulate (&DoFTools::template compute_intergrid_weights_3<dim>)
.collect_args (coarse_grid, coarse_component,
coarse_to_fine_grid_map, parameter_dofs,
weight_mapping, weights,
template <int dim>
void
-DoFTools::compute_intergrid_weights_1 (const DoFHandler<dim> &coarse_grid,
+DoFTools::compute_intergrid_weights_3 (const DoFHandler<dim> &coarse_grid,
const unsigned int coarse_component,
const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
const std::vector<Vector<double> > ¶meter_dofs,
const DoFHandler<deal_II_dimension> &,
const unsigned int ,
const InterGridMap<DoFHandler,deal_II_dimension> &,
- ConstraintMatrix &,
- std::vector<std::map<unsigned int, float> > *);
+ ConstraintMatrix &);
+template
+void
+DoFTools::compute_intergrid_transfer_representation (const DoFHandler<deal_II_dimension> &,
+ const unsigned int ,
+ const DoFHandler<deal_II_dimension> &,
+ const unsigned int ,
+ const InterGridMap<DoFHandler,deal_II_dimension> &,
+ std::vector<std::map<unsigned int, float> > &);
projects / dealii.git / commitdiff