--- /dev/null
+/*---------------------------- solutiontransfer.h ----------------------*/
+/* $Id$ */
+/* Ralf Hartmann, University of Heidelberg */
+#ifndef __solutiontransfer_H
+#define __solutiontransfer_H
+/*---------------------------- solutiontransfer.h ----------------------*/
+
+#include <base/exceptions.h>
+#include <vector.h>
+template <typename number> class Vector;
+template <int dim> class DoFHandler;
+
+/**
+ * Transfers a discrete FE function (like a solution vector) by interpolation
+ * while refining and/or coarsening a grid. During interpolation the
+ * vector is reinitialized to the new size and filled with the interpolated
+ * values.
+ *
+ * \subsection{Usage}
+ *
+ * As the interpolation while
+ * coarsening is much more complicated to organize
+ * (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
+ * (i.e. not locally coarsened) then use #SolutionTransfer# as follows
+ * \begin{verbatim}
+ * SolutionTransfer<dim, double> soltrans(*dof_handler);
+ * soltrans.prepare_for_pure_refinement();
+ * // some refinement e.g.
+ * tria->refine_and_coarsen_fixed_fraction(error_indicator, 0.3, 0);
+ * tria->execute_coarsening_and_refinement();
+ * soltrans.refine_interpolate(solution);
+ * // if necessary interpolate some
+ * // more functions
+ * soltrans.refine_interpolate(sol2);
+ * ...
+ * \end{verbatim}
+ * \item If the grid will be coarsenend and refined
+ * then use #SolutionTransfer# as follows
+ * \begin{verbatim}
+ * SolutionTransfer<dim, double> soltrans(*dof_handler);
+ * // some refinement e.g.
+ * tria->refine_and_coarsen_fixed_fraction(error_indicator, 0.3, 0);
+ * // very important:
+ * tria->prepare_coarsening_and_refinement();
+ * soltrans.prepare_for_coarsening_and_refinement(solution);
+ * tria->execute_coarsening_and_refinement ();
+ * soltrans.refine_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.
+ * \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
+ * them otherwise.
+ *
+ * @author Ralf Hartmann, 1999
+ */
+template<int dim, typename number>
+class SolutionTransfer
+{
+ public:
+ /**
+ * Constructor, takes the current #DoFHandler#
+ * as argument.
+ */
+ SolutionTransfer(const DoFHandler<dim> &dof);
+
+ /**
+ * Destructor
+ */
+ ~SolutionTransfer();
+
+ /**
+ * Reinit this class to the state that
+ * it has
+ * directly after calling the Constructor
+ */
+ void reinit();
+
+ /**
+ * Prepares the #SolutionTransfer# for
+ * pure refinement. It
+ * stores the dof indices of each cell.
+ * After calling this function
+ * only calling the #refine_interpolate#
+ * functions is allowed.
+ */
+ void prepare_for_pure_refinement();
+
+ /**
+ * Prepares the #SolutionTransfer# for
+ * coarsening and refinement. It
+ * stores the dof indices of each cell and
+ * stores the dof values of the vectors in
+ * #all_in# in each cell that'll be coarsened.
+ * #all_in# includes all vectors
+ * that are to be interpolated
+ * onto the new (refined and/or
+ * coarsenend) grid.
+ */
+ void prepare_for_coarsening_and_refinement (const vector<Vector<number> > &all_in);
+
+ /**
+ * Same as previous function
+ * but only for 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.
+ * It assumes the vectors having the
+ * right sizes (i.e. in.size()==n_dofs_old,
+ * out.size()==n_dofs_refined)
+ *
+ * Calling this function is allowed only
+ * if #prepare_for_pure_refinement# is called
+ * and the refinement is
+ * executed before.
+ * Multiple calling of this function is
+ * allowed. e.g. for interpolating several
+ * functions.
+ */
+ void refine_interpolate (const Vector<number> &in,
+ Vector<number> &out) const;
+
+ /**
+ * Same as #interpolate(Vector<number> in,
+ * Vector<number> out)# but it interpolates
+ * just 'in-place'. Therefore #vec# will be
+ * reinitialized to the new needed vector
+ * dimension.
+ */
+ void refine_interpolate (Vector<number> &vec) const;
+
+ /**
+ * This function
+ * 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
+ * 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).
+ *
+ * Calling this function is allowed only
+ * if first #Triangulation::prepare_coarsening_
+ * and_refinement#, second
+ * #SolutionTransfer::prepare_for_coarsening_
+ * and_refinement#, an then third
+ * #Triangulation::execute_coarsening_
+ * and_refinement# are called before.
+ * Multiple calling of this function is
+ * NOT allowed. Interpolating
+ * several functions can be performed
+ * in one step.
+ */
+ void interpolate (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
+ * several functions can be performed
+ * in one step by using #void
+ * interpolate (vector<Vector<number>
+ * >&all_out) const#
+ */
+ void interpolate (Vector<number> &out) const;
+
+ /**
+ * Exception
+ */
+ DeclException0(ExcNotPrepared);
+
+ /**
+ * Exception
+ */
+ DeclException0(ExcAlreadyPrepForRef);
+
+ /**
+ * Exception
+ */
+ DeclException0(ExcAlreadyPrepForRefAndCoarse);
+
+ /**
+ * Exception
+ */
+ DeclException0(ExcTriaPrepCoarseningNotCalledBefore);
+
+ /**
+ * Exception
+ */
+ DeclException0(ExcNoInVectorsGiven);
+
+ /**
+ * Exception
+ */
+ DeclException0(ExcVectorsDifferFromInVectors);
+
+ /**
+ * Exception
+ */
+ DeclException2(ExcWrongVectorSize,
+ int, int,
+ << "The size of the vector is " << arg1
+ << "although it should be " << arg2 << ".");
+
+ /**
+ * Exception
+ */
+ DeclException0(ExcInternalError);
+
+ private:
+
+ /**
+ * Pointer to the degree of freedom handler
+ * to work with.
+ */
+ const DoFHandler<dim> *dof_handler;
+
+ /**
+ * Stores the number of DoFs before the
+ * refinement and/or coarsening.
+ */
+ unsigned int n_dofs_old;
+
+ /**
+ * Denotes whether the #SolutionTransfer#
+ * is 'prepared for pure refinement'
+ * or not.
+ */
+ bool prepared_for_pure_refinement;
+
+ /**
+ * Denotes whether the #SolutionTransfer#
+ * is 'prepared for coarsening and refinement'
+ * or not.
+ */
+ bool prepared_for_coarsening_and_refinement;
+
+ /**
+ * Is used for #prepare_for_refining#
+ * (of course also for
+ * #repare_for_refining_and_coarsening#)
+ * and stores all dof indices
+ * of the cells that'll be refined
+ */
+ vector<vector<int> > indices_on_cell;
+
+ /**
+ * All cell data (the dof indices and
+ * the dof values)
+ * should be accessable from each cell.
+ * As each cell has got only one
+ * #user_pointer#, multiple pointers to the
+ * data need to be packetized in a structure.
+ * Note that in our case on each cell
+ * either the
+ * #vector<int> indices# (if the cell
+ * will be refined) or the
+ * #vector<double> dof_values# (if the
+ * children of this cell will be deleted)
+ * is needed, hence one user_pointer should
+ * be sufficient, but to allow some errorchecks
+ * and to preserve the user from making
+ * user errors the #user_pointer# will be
+ * 'multiplied' by this structure.
+ */
+ struct Pointerstruct {
+ vector<int> *indices_ptr;
+ vector<Vector<number> > *dof_values_ptr;
+ };
+
+ /**
+ * Vector of all #Pointerstructs# (cf. there).
+ * It makes it
+ * easier to delete all these structs
+ * (without going over all #cell->user_pointer#)
+ * after they are not used any more, and
+ * collecting all these structures in a vector
+ * helps avoiding fraqmentation of the memory.
+ */
+ vector<Pointerstruct> all_pointerstructs;
+
+ /**
+ * Is used for
+ * #prepare_for_refining_and_coarsening#
+ * The interpolated dof values
+ * of all cells that'll be coarsened
+ * will be stored in this vector.
+ */
+ vector<vector<Vector<number> > > dof_values_on_cell;
+
+ /**
+ * After calling
+ * #prepare_for_refinement_and_coarsening
+ * (vector<Vector<number> > &all_in)#
+ * this this pointer points to the vector
+ * #all_in# for later comparison with
+ * the vector #all_out#
+ */
+ const vector<Vector<number> > * vecs_ptr;
+};
+
+
+
+
+
+
+/*---------------------------- solutiontransfer.h ---------------------------*/
+/* end of #ifndef __solutiontransfer_H */
+#endif
+/*---------------------------- solutiontransfer.h ---------------------------*/
--- /dev/null
+/*---------------------------- solutiontransfer.cc ---------------------*/
+/* $Id$ */
+/* Ralf Hartmann, University of Heidelberg */
+
+#include <grid/tria.h>
+#include <grid/dof.h>
+#include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
+#include <fe/fe.h>
+#include <lac/vector.h>
+#include <numerics/solutiontransfer.h>
+
+
+
+template<int dim, typename number>
+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) {}
+
+
+template<int dim, typename number>
+SolutionTransfer<dim, number>::~SolutionTransfer()
+{
+ reinit();
+}
+
+
+template<int dim, typename number>
+void SolutionTransfer<dim, number>::reinit()
+{
+ if (indices_on_cell.size())
+ indices_on_cell.erase(indices_on_cell.begin(), indices_on_cell.end());
+ if (all_pointerstructs.size())
+ 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;
+}
+
+
+
+template<int dim, typename number>
+void SolutionTransfer<dim, number>::prepare_for_pure_refinement()
+{
+ Assert(!prepared_for_pure_refinement, ExcAlreadyPrepForRef());
+ Assert(!prepared_for_coarsening_and_refinement,
+ ExcAlreadyPrepForRefAndCoarse());
+
+ reinit();
+
+ unsigned int n_cells=dof_handler->get_tria().n_active_cells();
+ 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));
+
+ DoFHandler<dim>::cell_iterator cell = dof_handler->begin(),
+ endc = dof_handler->end();
+
+ for (unsigned int i=0; cell!=endc; ++cell)
+ {
+ if (cell->active())
+ {
+ cell->set_user_pointer(&indices_on_cell[i]);
+ cell->get_dof_indices(indices_on_cell[i]);
+ ++i;
+ }
+ else
+ cell->clear_user_pointer();
+ }
+ prepared_for_pure_refinement=true;
+}
+
+
+template<int dim, typename number>
+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));
+ Assert(out.size()==dof_handler->n_dofs(),
+ ExcWrongVectorSize(out.size(),dof_handler->n_dofs()));
+
+ unsigned int total_dofs=dof_handler->get_fe().total_dofs;
+ Vector<number> local_values(total_dofs);
+
+ DoFHandler<dim>::cell_iterator cell = dof_handler->begin(),
+ endc = dof_handler->end();
+
+ vector<int> *indexptr;
+
+ for (; cell!=endc; ++cell)
+ {
+ if (cell->user_pointer())
+ {
+ indexptr=static_cast<vector<int> *>(cell->user_pointer());
+ for (unsigned int i=0; i<total_dofs; ++i)
+ local_values(i)=in(indexptr->operator[](i));
+ cell->set_dof_values_by_interpolation(local_values, out);
+ }
+ }
+}
+
+template<int dim, typename number>
+void SolutionTransfer<dim, number>::refine_interpolate(
+ Vector<number> &vec) const
+{
+ Assert(vec.size()==n_dofs_old, ExcWrongVectorSize(vec.size(),n_dofs_old));
+
+ Vector<number> vec_old(vec);
+ vec.reinit(dof_handler->n_dofs());
+
+ refine_interpolate(vec_old, vec);
+}
+
+
+
+template<int dim, typename number>
+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();
+ Assert(in_size!=0, ExcNoInVectorsGiven());
+
+ reinit();
+
+ unsigned int n_cells=dof_handler->get_tria().n_active_cells();
+ 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;
+
+ // sollte vorher passiert sein
+// dof_handler->get_tria().prepare_coarsening();
+
+
+ // first count the number
+ // of cells that'll be coarsened
+ // and that'll stay or be refined
+ unsigned int n_cells_to_coarsen=0;
+ unsigned int n_cells_to_stay_or_refine=0;
+ DoFHandler<dim>::active_cell_iterator act_cell = dof_handler->begin_active(),
+ endc = dof_handler->end();
+ for (; act_cell!=endc; ++act_cell)
+ {
+ if (act_cell->coarsen_flag_set())
+ ++n_cells_to_coarsen;
+ else
+ ++n_cells_to_stay_or_refine;
+ }
+ Assert((n_cells_to_coarsen+n_cells_to_stay_or_refine)==n_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;
+
+ // 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)));
+// dof_values_on_cell=<vector<Vector<number> > (in_size, Vector<number> (total_dofs))> (vector<vector<Vector<number> > >, n_coarsen_fathers);
+
+
+ cout << "dof_values_on_cell.size()=" << dof_values_on_cell.size()
+ << ", n_coarsen_fathers=" << n_coarsen_fathers << endl;
+ if (dof_values_on_cell.size())
+ {
+ cout << endl << "dof_values_on_cell[0].size()=" << dof_values_on_cell[0].size()
+ << ", in_size=" << in_size << endl;
+ cout << "dof_values_on_cell[0][0].size()=" << dof_values_on_cell[0][0].size()
+ << ", total_dofs=" << total_dofs << endl;
+ }
+
+ all_pointerstructs=vector<Pointerstruct> (
+ n_cells_to_stay_or_refine+n_coarsen_fathers);
+
+
+ // we need counters for
+ // the 'to_stay_or_refine' cells 'n_sr',
+ // the 'coarsen_fathers' cells 'n_cf',
+ // and all the cells where a
+ // #Pointerstruct# is needed 'n'
+ unsigned int n_sr=0, n_cf=0, n=0;
+// Vector<number> values(total_dofs);
+ DoFHandler<dim>::cell_iterator cell = dof_handler->begin();
+ for (; cell!=endc; ++cell)
+ {
+ if (cell->active() && !cell->coarsen_flag_set())
+ {
+ 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;
+ }
+ else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
+ {
+ 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]);
+ }
+ 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;
+ }
+ else
+ cell->clear_user_pointer();
+ }
+ Assert(n_sr==n_cells_to_stay_or_refine, ExcInternalError());
+ Assert(n_cf==n_coarsen_fathers, 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)
+{
+ 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
+{
+ Assert(prepared_for_coarsening_and_refinement, ExcNotPrepared());
+ Assert(vecs_ptr==&all_out, ExcVectorsDifferFromInVectors());
+ 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;
+ Vector<number> local_values(total_dofs);
+
+ vector<int> *indexptr;
+ Pointerstruct *structptr;
+ vector<Vector<number> > *valuesptr;
+ vector<int> dofs(total_dofs);
+
+ DoFHandler<dim>::cell_iterator cell = dof_handler->begin(),
+ endc = dof_handler->end();
+ unsigned int counter=0;
+
+ for (; cell!=endc; ++cell)
+ {
+ if (cell->user_pointer())
+ {
+ structptr=static_cast<Pointerstruct *>(cell->user_pointer());
+ // cell stayed or is refined
+ if (indexptr=structptr->indices_ptr)
+ {
+ 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]);
+ }
+ }
+ // children of cell were deleted
+ else if (valuesptr=structptr->dof_values_ptr)
+ {
+// Assert(!cell->has_children(), ExcInternalError());
+ if (cell->has_children())
+ cout << "counter=" << ++counter << endl;
+ cell->get_dof_indices(dofs);
+ for (unsigned int j=0; j<out_size; ++j)
+ {
+ for (unsigned int i=0; i<total_dofs; ++i)
+ all_out[j](dofs[i])=valuesptr->operator[](j)(i);
+ }
+ }
+ // undefined status
+ else
+ Assert(false, ExcInternalError());
+ }
+ }
+ // delete all stored data
+ // should be called outside
+// reinit();
+}
+
+
+template<int dim, typename number>
+void SolutionTransfer<dim, number>::interpolate(Vector<number> &out) const
+{
+ vector<Vector<number> > all_in=vector<Vector<number> >(1, out);
+ interpolate(all_in);
+ out=all_in[0];
+}
+
+
+template class SolutionTransfer<2, double>;
+
+/*---------------------------- solutiontransfer.cc ----------------------*/