//---------------------------- dof_tools.cc ---------------------------
+#include <base/multithread_info.h>
+#include <base/thread_management.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <grid/intergrid_map.h>
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
};
+ // 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
+ //
+ // do this in a separate function
+ // to allow for multithreading
+ // there. see this function also if
+ // 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);
+
+
+ // 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. the
+ // representant will become an
+ // unconstrained degree of freedom,
+ // while all others will be
+ // constrained to this dof (and
+ // possibly others)
+ 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;
+ };
+
+
+ // 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...
+ vector<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
+ {
+ 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);
+
+ constraint_line.clear ();
+ for (unsigned int row=first_used_row; row<weights.m(); ++row)
+ if (weights(row,weight_mapping[global_dof]) != 0)
+ constraint_line.push_back (make_pair(representants[row],
+ weights(row,weight_mapping[global_dof])));
+
+ constraints.add_entries (global_dof, constraint_line);
+ };
+};
+
+
+
+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 vector<Vector<double> > ¶meter_dofs,
+ const vector<int> &weight_mapping,
+ FullMatrix<double> &weights)
+{
+ // simply distribute the range of
+ // cells to different threads
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ vector<pair<active_cell_iterator,active_cell_iterator> >
+ cell_intervals = Threads::split_range<active_cell_iterator> (coarse_grid.begin_active(),
+ coarse_grid.end(),
+ multithread_info.n_default_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>)
+ .collect_args (coarse_grid, coarse_component,
+ coarse_to_fine_grid_map, parameter_dofs,
+ weight_mapping, weights,
+ cell_intervals[i].first,
+ cell_intervals[i].second));
+
+ // wait for the threads to finish
+ thread_manager.wait ();
+};
+
+
+
+template <int dim>
+void
+DoFTools::compute_intergrid_weights_1 (const DoFHandler<dim> &coarse_grid,
+ const unsigned int coarse_component,
+ const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+ const vector<Vector<double> > ¶meter_dofs,
+ const vector<int> &weight_mapping,
+ FullMatrix<double> &weights,
+ const typename DoFHandler<dim>::active_cell_iterator &begin,
+ const typename DoFHandler<dim>::active_cell_iterator &end)
+{
+ // aliases to the finite elements
+ // used by the dof handlers:
+ const FiniteElement<dim> &coarse_fe = coarse_grid.get_fe();
+
// for each cell on the parameter grid:
// find out which degrees of freedom on the
// fine grid correspond in which way to
// both cells 1 and 2, but the
// correct weight is nevertheless
// only 1.
+
+ // vector to hold the representation of
+ // a single degree of freedom on the
+ // coarse grid (for the selected fe)
+ // on the fine grid
+ const unsigned int n_fine_dofs = weight_mapping.size();
+ Vector<double> global_parameter_representation (n_fine_dofs);
+
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)
+ for (cell=begin; cell!=end; ++cell)
{
// get the global indices of the
// parameter dofs on this parameter
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. the
- // representant will become an
- // unconstrained degree of freedom,
- // while all others will be
- // constrained to this dof (and
- // possibly others)
- 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]));
- };
};
void
DoFTools::distribute_cell_to_dof_vector (const DoFHandler<deal_II_dimension> &dof_handler,
const Vector<float> &cell_data,
- Vector<double> &dof_data) const;
+ Vector<double> &dof_data);
template
void
DoFTools::distribute_cell_to_dof_vector (const DoFHandler<deal_II_dimension> &dof_handler,
const Vector<double> &cell_data,
- Vector<double> &dof_data) const;
+ Vector<double> &dof_data);
template void DoFTools::extract_dofs(const DoFHandler<deal_II_dimension>& dof,
projects / dealii-svn.git / commitdiff