*
* As the interpolation while
* coarsening is much more complicated to organize
- * (see further documentation below) than interpolation while pure refinement
+ * (see further documentation below) than interpolation while pure refinement,
* #SolutionTransfer# offers two possible usages.
* \begin{itemize}
* \item If the grid will only be purely refined
* // some refinement e.g.
* tria->refine_and_coarsen_fixed_fraction(error_indicator, 0.3, 0);
* tria->execute_coarsening_and_refinement();
+ * dof_handler->distribute_dofs (fe);
* soltrans.refine_interpolate(solution);
* // if necessary interpolate some
* // more functions
* tria->prepare_coarsening_and_refinement();
* soltrans.prepare_for_coarsening_and_refinement(solution);
* tria->execute_coarsening_and_refinement ();
- * soltrans.refine_interpolate(solution);
+ * dof_handler->distribute_dofs (fe);
+ * soltrans.interpolate(solution);
* \end{verbatim}
+
* Multiple calling of this function #void interpolate (Vector<number> &out) const#
- * is NOT allowed. Interpolating
- * several functions can be performed
- * in one step by using #void
- * interpolate (vector<Vector<number>
- * >&all_out) const#, see there.
+ * is NOT allowed. Interpolating several functions can be performed in one step
+ * by using #void interpolate (vector<Vector<number> >&all_out) const#, and
+ * using the respective #prepare_for_coarsening_and_refinement# function taking
+ * several vectors as input before actually refining and coarsening the
+ * triangulation (see there).
* \end{itemize}
- * For deleting all stored data in and reinitializing the #SolutionTransfer#
- * use the #reinit()# function.
*
- * Note that the #user_pointer# of the cell are used. Be sure that you don't need
+ * For deleting all stored data in and reinitializing the
+ * #SolutionTransfer# use the #clear()# function.
+ *
+ * Note that the #user_pointer# of some cells are used. Be sure that you don't need
* them otherwise.
*
- * @author Ralf Hartmann, 1999
- */
+ * The template argument #number# denotes the data type of the vectors you want
+ * to transfer.
+ *
+ * @author Ralf Hartmann, 1999 */
template<int dim, typename number>
class SolutionTransfer
{
* it has
* directly after calling the Constructor
*/
- void reinit();
+ void clear();
/**
* Prepares the #SolutionTransfer# for
/**
* Same as previous function
- * but only for one discrete function
+ * but for only one discrete function
* to interpolate.
*/
void prepare_for_coarsening_and_refinement (const Vector<number> &in);
/**
* This function
- * interpolates the discrete function #in#
- * on the grid before the
- * refinement to the function #out#
- * on the refined grid.
+ * interpolates the discrete function #in#,
+ * which is a vector on the grid before the
+ * refinement, to the function #out#
+ * which then is a vector on the refined grid.
* It assumes the vectors having the
* right sizes (i.e. in.size()==n_dofs_old,
* out.size()==n_dofs_refined)
* interpolates the discrete functions
* that are stored in #all_out# onto
* the refined and/or coarsenend grid.
- * It assumes the vectors in #all_out#
- * being idendically the same vectors
+ * It assumes the vectors in #all_in#
+ * denote the same vectors
* as in #all_in# as parameter
* of #prepare_for_refinement_and_coarsening
* (vector<Vector<number> > &all_in)#.
- * So they have the
- * right sizes (i.e. all_out[j].size()==
- * n_dofs_old for all j).
+ * However, there is no way of verifying
+ * this internally, so be careful here.
*
* Calling this function is allowed only
* if first #Triangulation::prepare_coarsening_
* several functions can be performed
* in one step.
*/
- void interpolate (vector<Vector<number> >&all_out) const;
+ void interpolate (const vector<Vector<number> >&all_in,
+ vector<Vector<number> > &all_out) const;
/**
* Same as the previous function.
* It interpolates only one function.
- * It assumes the vector having the
- * right size (i.e. all_out[j].size()==
- * n_dofs_old for all j)
*
* Multiple calling of this function is
* NOT allowed. Interpolating
* interpolate (vector<Vector<number>
* >&all_out) const#
*/
- void interpolate (Vector<number> &out) const;
+ void interpolate (const Vector<number> &in,
+ Vector<number> &out) const;
/**
* Exception
* Exception
*/
DeclException0(ExcInternalError);
-
+
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidVectorSize,
+ int, int,
+ << "The data vector has the size " << arg1
+ << ", but " << arg2 << " was expected.");
+
private:
/**
* After calling
* #prepare_for_refinement_and_coarsening
* (vector<Vector<number> > &all_in)#
- * this this pointer points to the vector
+ * this pointer points to the vector
* #all_in# for later comparison with
* the vector #all_out#
*/
- const vector<Vector<number> > * vecs_ptr;
+// const vector<Vector<number> > * vecs_ptr;
};
template<int dim, typename number>
-SolutionTransfer<dim, number>::SolutionTransfer(
- const DoFHandler<dim> &dof):
+SolutionTransfer<dim, number>::SolutionTransfer(const DoFHandler<dim> &dof):
dof_handler(&dof),
n_dofs_old(0),
prepared_for_pure_refinement(0),
- prepared_for_coarsening_and_refinement(0),
- vecs_ptr(0) {}
+ prepared_for_coarsening_and_refinement(0)
+{}
+
template<int dim, typename number>
SolutionTransfer<dim, number>::~SolutionTransfer()
{
- reinit();
+ clear ();
}
+
template<int dim, typename number>
-void SolutionTransfer<dim, number>::reinit()
+void SolutionTransfer<dim, number>::clear ()
{
if (indices_on_cell.size())
indices_on_cell.erase(indices_on_cell.begin(), indices_on_cell.end());
all_pointerstructs.erase(all_pointerstructs.begin(), all_pointerstructs.end());
if (dof_values_on_cell.size())
dof_values_on_cell.erase(dof_values_on_cell.begin(), dof_values_on_cell.end());
- vecs_ptr=0;
+
prepared_for_pure_refinement=false;
prepared_for_coarsening_and_refinement=false;
}
Assert(!prepared_for_coarsening_and_refinement,
ExcAlreadyPrepForRefAndCoarse());
- reinit();
+ clear();
- unsigned int n_cells=dof_handler->get_tria().n_active_cells();
- unsigned int total_dofs=dof_handler->get_fe().total_dofs;
+ const unsigned int n_active_cells = dof_handler->get_tria().n_active_cells();
+ const unsigned int total_dofs = dof_handler->get_fe().total_dofs;
n_dofs_old=dof_handler->n_dofs();
- indices_on_cell=vector<vector<int> > (
- n_cells, vector<int> (total_dofs));
+ indices_on_cell=vector<vector<int> > (n_active_cells,
+ vector<int> (total_dofs));
DoFHandler<dim>::cell_iterator cell = dof_handler->begin(),
endc = dof_handler->end();
{
if (cell->active())
{
- cell->set_user_pointer(&indices_on_cell[i]);
+ // on each cell store the indices
+ // of the dofs. after refining we
+ // get the values on the children
+ // by taking these indices, getting
+ // the respective values out of
+ // the data vectors and prolonging
+ // them to the children
cell->get_dof_indices(indices_on_cell[i]);
+ cell->set_user_pointer(&indices_on_cell[i]);
+
++i;
}
else
}
+
template<int dim, typename number>
-void SolutionTransfer<dim, number>::refine_interpolate(
- const Vector<number> &in, Vector<number> &out) const
+void
+SolutionTransfer<dim, number>::refine_interpolate(const Vector<number> &in,
+ Vector<number> &out) const
{
Assert(prepared_for_pure_refinement, ExcNotPrepared());
Assert(in.size()==n_dofs_old, ExcWrongVectorSize(in.size(),n_dofs_old));
for (; cell!=endc; ++cell)
{
if (cell->user_pointer())
+ // this cell was refined or not
+ // touched at all, so we can get
+ // the new values by just setting
+ // or interpolating to the children,
+ // which is both done by one
+ // function
{
indexptr=static_cast<vector<int> *>(cell->user_pointer());
for (unsigned int i=0; i<total_dofs; ++i)
}
}
+
+
template<int dim, typename number>
-void SolutionTransfer<dim, number>::refine_interpolate(
- Vector<number> &vec) const
+void SolutionTransfer<dim, number>::refine_interpolate (Vector<number> &vec) const
{
Assert(vec.size()==n_dofs_old, ExcWrongVectorSize(vec.size(),n_dofs_old));
template<int dim, typename number>
-void SolutionTransfer<dim, number>::prepare_for_coarsening_and_refinement(
- const vector<Vector<number> > &all_in)
+void
+SolutionTransfer<dim, number>::
+prepare_for_coarsening_and_refinement(const vector<Vector<number> > &all_in)
{
Assert(!prepared_for_pure_refinement, ExcAlreadyPrepForRef());
Assert(!prepared_for_coarsening_and_refinement,
ExcAlreadyPrepForRefAndCoarse());
- unsigned int in_size=all_in.size();
+
+ const unsigned int in_size=all_in.size();
Assert(in_size!=0, ExcNoInVectorsGiven());
- reinit();
+ clear();
- unsigned int n_cells=dof_handler->get_tria().n_active_cells();
- unsigned int total_dofs=dof_handler->get_fe().total_dofs;
+ const unsigned int n_active_cells = dof_handler->get_tria().n_active_cells();
+ const unsigned int total_dofs = dof_handler->get_fe().total_dofs;
n_dofs_old=dof_handler->n_dofs();
for (unsigned int i=0; i<in_size; ++i)
Assert(all_in[i].size()==n_dofs_old,
ExcWrongVectorSize(all_in[i].size(),n_dofs_old));
}
- vecs_ptr=&all_in;
+
// first count the number
// of cells that'll be coarsened
// and that'll stay or be refined
else
++n_cells_to_stay_or_refine;
}
- Assert((n_cells_to_coarsen+n_cells_to_stay_or_refine)==n_cells,
+ Assert((n_cells_to_coarsen+n_cells_to_stay_or_refine)==n_active_cells,
ExcInternalError());
Assert(n_cells_to_coarsen%GeometryInfo<dim>::children_per_cell==0,
ExcInternalError());
Assert(n_cells_to_coarsen%GeometryInfo<dim>::children_per_cell==0,
ExcTriaPrepCoarseningNotCalledBefore());
- unsigned int n_coarsen_fathers=n_cells_to_coarsen/
- GeometryInfo<dim>::children_per_cell;
+ const unsigned int n_coarsen_fathers = n_cells_to_coarsen /
+ GeometryInfo<dim>::children_per_cell;
// allocate the needed memory
- indices_on_cell=vector<vector<int> > (
- n_cells_to_stay_or_refine, vector<int> (total_dofs));
- dof_values_on_cell=vector<vector<Vector<number> > > (
- n_coarsen_fathers, vector<Vector<number> > (
- in_size, Vector<number> (total_dofs)));
+ indices_on_cell = vector<vector<int> > (n_cells_to_stay_or_refine,
+ vector<int> (total_dofs));
+ dof_values_on_cell
+ = vector<vector<Vector<number> > > (n_coarsen_fathers,
+ vector<Vector<number> > (in_size,
+ Vector<number> (total_dofs)));
- all_pointerstructs=vector<Pointerstruct> (
- n_cells_to_stay_or_refine+n_coarsen_fathers);
+ all_pointerstructs = vector<Pointerstruct> (n_cells_to_stay_or_refine +
+ n_coarsen_fathers);
// we need counters for
{
if (cell->active() && !cell->coarsen_flag_set())
{
+ // cell will not be coarsened,
+ // so we get away by storing the
+ // dof indices and later
+ // interpolating to the children
cell->get_dof_indices(indices_on_cell[n_sr]);
all_pointerstructs[n].indices_ptr=&indices_on_cell[n_sr];
all_pointerstructs[n].dof_values_ptr=0;
cell->set_user_pointer(&all_pointerstructs[n]);
- ++n_sr, ++n;
+ ++n_sr;
+ ++n;
}
else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
{
+ // note that if one child has the
+ // coaresen flag, then all should
+ // have if Tria::prepare_* has
+ // worked correctly
for (unsigned int i=1; i<GeometryInfo<dim>::children_per_cell; ++i)
Assert(cell->child(i)->coarsen_flag_set(),
ExcTriaPrepCoarseningNotCalledBefore());
for (unsigned int j=0; j<in_size; ++j)
{
- cell->get_interpolated_dof_values(
- all_in[j], dof_values_on_cell[n_cf][j]);
+ // store the data of each of
+ // the input vectors
+ cell->get_interpolated_dof_values(all_in[j],
+ dof_values_on_cell[n_cf][j]);
}
all_pointerstructs[n].indices_ptr=0;
all_pointerstructs[n].dof_values_ptr=&dof_values_on_cell[n_cf];
cell->set_user_pointer(&all_pointerstructs[n]);
- ++n_cf, ++n;
+ ++n_cf;
+ ++n;
}
else
+ // some cell on the lower levels to
+ // which nothing will happen
cell->clear_user_pointer();
}
Assert(n_sr==n_cells_to_stay_or_refine, ExcInternalError());
prepared_for_coarsening_and_refinement=true;
}
+
+
template<int dim, typename number>
-void SolutionTransfer<dim, number>::prepare_for_coarsening_and_refinement(
- const Vector<number> &in)
+void
+SolutionTransfer<dim, number>::prepare_for_coarsening_and_refinement(const Vector<number> &in)
{
vector<Vector<number> > all_in=vector<Vector<number> >(1, in);
prepare_for_coarsening_and_refinement(all_in);
template<int dim, typename number>
-void SolutionTransfer<dim, number>::interpolate (
- vector<Vector<number> > &all_out) const
+void SolutionTransfer<dim, number>::
+interpolate (const vector<Vector<number> > &all_in,
+ vector<Vector<number> > &all_out) const
{
Assert(prepared_for_coarsening_and_refinement, ExcNotPrepared());
- Assert(vecs_ptr==&all_out, ExcVectorsDifferFromInVectors());
+ for (unsigned int i=0; i<all_in.size(); ++i)
+ Assert (all_in[i].size() == n_dofs_old,
+ ExcInvalidVectorSize(all_in[i].size(), n_dofs_old));
+
unsigned int out_size=all_out.size();
- // local copy of the vectors
- vector<Vector<number> > all_in(all_out);
- // resize the out_all vectors
- for (unsigned int j=0; j<out_size; ++j)
- all_out[j].reinit(dof_handler->n_dofs());
- unsigned int total_dofs=dof_handler->get_fe().total_dofs;
+ // resize the output vector
+ if (all_out.size() != all_in.size())
+ all_out = vector<Vector<number> >(all_in.size(),
+ Vector<number>(dof_handler->n_dofs()));
+ else
+ {
+ for (unsigned int i=0; i<all_in.size(); ++i)
+ if (all_out[i].size() != all_in[i].size())
+ all_out[i].reinit (all_in[i].size());
+ };
+
+
+ const unsigned int total_dofs=dof_handler->get_fe().total_dofs;
Vector<number> local_values(total_dofs);
- vector<int> *indexptr;
- Pointerstruct *structptr;
+ vector<int> *indexptr;
+ Pointerstruct *structptr;
vector<Vector<number> > *valuesptr;
- vector<int> dofs(total_dofs);
+ vector<int> dofs(total_dofs);
DoFHandler<dim>::cell_iterator cell = dof_handler->begin(),
endc = dof_handler->end();
// cell stayed or is refined
if (indexptr=structptr->indices_ptr)
{
+ Assert (structptr->dof_values_ptr == 0,
+ ExcInternalError());
+
+ // get the values of each of
+ // the input data vectors on
+ // this cell and prolong it
+ // to its children
for (unsigned int j=0; j<out_size; ++j)
{
for (unsigned int i=0; i<total_dofs; ++i)
local_values(i)=all_in[j](indexptr->operator[](i));
- cell->set_dof_values_by_interpolation(
- local_values, all_out[j]);
+ cell->set_dof_values_by_interpolation(local_values,
+ all_out[j]);
}
}
// children of cell were deleted
else if (valuesptr=structptr->dof_values_ptr)
{
- Assert(!cell->has_children(), ExcInternalError());
+ Assert (!cell->has_children(), ExcInternalError());
+ Assert (structptr->indices_ptr == 0,
+ ExcInternalError());
+
+ // get the local indices
cell->get_dof_indices(dofs);
+
+ // distribute the stored data
+ // to the new vectors
for (unsigned int j=0; j<out_size; ++j)
{
for (unsigned int i=0; i<total_dofs; ++i)
}
+
template<int dim, typename number>
-void SolutionTransfer<dim, number>::interpolate(Vector<number> &out) const
+void SolutionTransfer<dim, number>::interpolate(const Vector<number> &in,
+ Vector<number> &out) const
{
- vector<Vector<number> > all_in=vector<Vector<number> >(1, out);
- interpolate(all_in);
- out=all_in[0];
+ vector<Vector<number> > all_in=vector<Vector<number> >(1, in);
+ vector<Vector<number> > all_out=vector<Vector<number> >(1, out);
+ interpolate(all_in,
+ all_out);
+ out=all_out[0];
}
-template class SolutionTransfer<2, double>;
+
+
+template class SolutionTransfer<deal_II_dimension, float>;
+template class SolutionTransfer<deal_II_dimension, double>;
/*---------------------------- solutiontransfer.cc ----------------------*/