From 8d28b7cd3669612bb31c5aea521b7a1bea01735f Mon Sep 17 00:00:00 2001 From: hartmann Date: Tue, 9 Mar 1999 15:41:18 +0000 Subject: [PATCH] Initial version of the SolutionTransfer class that interpolates discrete functions while refinement and/or coarsening. git-svn-id: https://svn.dealii.org/trunk@974 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/solution_transfer.h | 336 ++++++++++++++++++ .../source/numerics/solution_transfer.cc | 329 +++++++++++++++++ 2 files changed, 665 insertions(+) create mode 100644 deal.II/deal.II/include/numerics/solution_transfer.h create mode 100644 deal.II/deal.II/source/numerics/solution_transfer.cc diff --git a/deal.II/deal.II/include/numerics/solution_transfer.h b/deal.II/deal.II/include/numerics/solution_transfer.h new file mode 100644 index 0000000000..cdc974f23e --- /dev/null +++ b/deal.II/deal.II/include/numerics/solution_transfer.h @@ -0,0 +1,336 @@ +/*---------------------------- solutiontransfer.h ----------------------*/ +/* $Id$ */ +/* Ralf Hartmann, University of Heidelberg */ +#ifndef __solutiontransfer_H +#define __solutiontransfer_H +/*---------------------------- solutiontransfer.h ----------------------*/ + +#include +#include +template class Vector; +template 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 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 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 &out) const# + * is NOT allowed. Interpolating + * several functions can be performed + * in one step by using #void + * interpolate (vector + * >&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 +class SolutionTransfer +{ + public: + /** + * Constructor, takes the current #DoFHandler# + * as argument. + */ + SolutionTransfer(const DoFHandler &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 > &all_in); + + /** + * Same as previous function + * but only for one discrete function + * to interpolate. + */ + void prepare_for_coarsening_and_refinement (const Vector &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 &in, + Vector &out) const; + + /** + * Same as #interpolate(Vector in, + * Vector out)# but it interpolates + * just 'in-place'. Therefore #vec# will be + * reinitialized to the new needed vector + * dimension. + */ + void refine_interpolate (Vector &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 > &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 >&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 + * >&all_out) const# + */ + void interpolate (Vector &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 *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 > 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 indices# (if the cell + * will be refined) or the + * #vector 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 *indices_ptr; + vector > *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 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 > > dof_values_on_cell; + + /** + * After calling + * #prepare_for_refinement_and_coarsening + * (vector > &all_in)# + * this this pointer points to the vector + * #all_in# for later comparison with + * the vector #all_out# + */ + const vector > * vecs_ptr; +}; + + + + + + +/*---------------------------- solutiontransfer.h ---------------------------*/ +/* end of #ifndef __solutiontransfer_H */ +#endif +/*---------------------------- solutiontransfer.h ---------------------------*/ diff --git a/deal.II/deal.II/source/numerics/solution_transfer.cc b/deal.II/deal.II/source/numerics/solution_transfer.cc new file mode 100644 index 0000000000..061e5e24d9 --- /dev/null +++ b/deal.II/deal.II/source/numerics/solution_transfer.cc @@ -0,0 +1,329 @@ +/*---------------------------- solutiontransfer.cc ---------------------*/ +/* $Id$ */ +/* Ralf Hartmann, University of Heidelberg */ + +#include +#include +#include +#include +#include +#include +#include +#include + + + +template +SolutionTransfer::SolutionTransfer( + const DoFHandler &dof): + dof_handler(&dof), + n_dofs_old(0), + prepared_for_pure_refinement(0), + prepared_for_coarsening_and_refinement(0), + vecs_ptr(0) {} + + +template +SolutionTransfer::~SolutionTransfer() +{ + reinit(); +} + + +template +void SolutionTransfer::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 +void SolutionTransfer::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 > ( + n_cells, vector (total_dofs)); + + DoFHandler::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 +void SolutionTransfer::refine_interpolate( + const Vector &in, Vector &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 local_values(total_dofs); + + DoFHandler::cell_iterator cell = dof_handler->begin(), + endc = dof_handler->end(); + + vector *indexptr; + + for (; cell!=endc; ++cell) + { + if (cell->user_pointer()) + { + indexptr=static_cast *>(cell->user_pointer()); + for (unsigned int i=0; ioperator[](i)); + cell->set_dof_values_by_interpolation(local_values, out); + } + } +} + +template +void SolutionTransfer::refine_interpolate( + Vector &vec) const +{ + Assert(vec.size()==n_dofs_old, ExcWrongVectorSize(vec.size(),n_dofs_old)); + + Vector vec_old(vec); + vec.reinit(dof_handler->n_dofs()); + + refine_interpolate(vec_old, vec); +} + + + +template +void SolutionTransfer::prepare_for_coarsening_and_refinement( + const vector > &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; iget_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::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::children_per_cell==0, + ExcInternalError()); + Assert(n_cells_to_coarsen%GeometryInfo::children_per_cell==0, + ExcTriaPrepCoarseningNotCalledBefore()); + + unsigned int n_coarsen_fathers=n_cells_to_coarsen/ + GeometryInfo::children_per_cell; + + // allocate the needed memory + indices_on_cell=vector > ( + n_cells_to_stay_or_refine, vector (total_dofs)); + dof_values_on_cell=vector > > ( + n_coarsen_fathers, vector > ( + in_size, Vector (total_dofs))); +// dof_values_on_cell= > (in_size, Vector (total_dofs))> (vector > >, 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 ( + 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 values(total_dofs); + DoFHandler::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::children_per_cell; ++i) + Assert(cell->child(i)->coarsen_flag_set(), + ExcTriaPrepCoarseningNotCalledBefore()); + + for (unsigned int j=0; jget_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 +void SolutionTransfer::prepare_for_coarsening_and_refinement( + const Vector &in) +{ + vector > all_in=vector >(1, in); + prepare_for_coarsening_and_refinement(all_in); +} + + + +template +void SolutionTransfer::interpolate ( + vector > &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 > all_in(all_out); + // resize the out_all vectors + for (unsigned int j=0; jn_dofs()); + + unsigned int total_dofs=dof_handler->get_fe().total_dofs; + Vector local_values(total_dofs); + + vector *indexptr; + Pointerstruct *structptr; + vector > *valuesptr; + vector dofs(total_dofs); + + DoFHandler::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(cell->user_pointer()); + // cell stayed or is refined + if (indexptr=structptr->indices_ptr) + { + for (unsigned int j=0; joperator[](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; joperator[](j)(i); + } + } + // undefined status + else + Assert(false, ExcInternalError()); + } + } + // delete all stored data + // should be called outside +// reinit(); +} + + +template +void SolutionTransfer::interpolate(Vector &out) const +{ + vector > all_in=vector >(1, out); + interpolate(all_in); + out=all_in[0]; +} + + +template class SolutionTransfer<2, double>; + +/*---------------------------- solutiontransfer.cc ----------------------*/ -- 2.39.5