]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
New: parallel::distributed::Triangulation::save()/load() to store
authorheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 12 Oct 2011 18:52:57 +0000 (18:52 +0000)
committerheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 12 Oct 2011 18:52:57 +0000 (18:52 +0000)
the refinement information to disk. Also supports saving solution vectors
using the SolutionTransfer class.

git-svn-id: https://svn.dealii.org/trunk@24595 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/distributed/solution_transfer.h
deal.II/include/deal.II/distributed/tria.h
deal.II/source/distributed/solution_transfer.cc
deal.II/source/distributed/tria.cc

index 5fcd80ed7bc2d7416a27d5cebdc66a6070491790..9d9bf74773115e7392361dd21840ba694340a0ec 100644 (file)
@@ -34,7 +34,9 @@ inconvenience this causes.
 <h3>General</h3>
 
 <ol>
-<li> None so far.
+<li> New: parallel::distributed::Triangulation::save()/load() to store
+the refinement information to disk. Also supports saving solution vectors
+using the SolutionTransfer class. (Timo Heister, 2011/10/12)
 </ol>
 
 
index 2d1229813f7124a5f8a62a5e36ed291daf0ba881..cf30c9ed44f66e61d458b322efde6e3e622884a0 100644 (file)
@@ -32,30 +32,69 @@ namespace parallel
  * coarsening a distributed grid and
  * handles the necessary communication.
  *
- * <h3>Usage</h3>
+ * @note: It is important to note, that if you use more than one
+ * SolutionTransfer object at the same time, that the calls to
+ * prepare_*() and interpolate()/deserialize() need to be in the same order.
+ *
+ * <h3>Note on ghost elements</h3>
+ * In a parallel computation PETSc or Trilinos vector may contain ghost
+ * elements or not. For reading in information with
+ * prepare_for_coarsening_and_refinement() or prepare_serialization() you need
+ * to supply vectors with ghost elements, so that all locally_active elements
+ * can be read. On the other hand, ghosted vectors are generally not writable,
+ * so for calls to interpolate() or deserialize() you need to supply
+ * distributed vectors without ghost elements.
+ * 
+ * <h3>Transfering a solution</h3>
+ * Here VECTOR is your favorite vector type, e.g. PETScWrappers::MPI::Vector,
+ * TrilinosWrappers::MPI::Vector, or corresponding blockvectors.
  * @verbatim
- * SolutionTransfer<dim, Vector<double> > soltrans(dof_handler);
- *                                     // flag some cells for refinement
- *                                     // and coarsening, e.g.
- * GridRefinement::refine_and_coarsen_fixed_fraction(
- *   *tria, error_indicators, 0.3, 0.05);
- *                                     // prepare the triangulation,
- * tria->prepare_coarsening_and_refinement();
- *                                     // prepare the SolutionTransfer object
- *                                     // for coarsening and refinement and give
- *                                     // the solution vector that we intend to
- *                                     // interpolate later,
- * soltrans.prepare_for_coarsening_and_refinement(solution);
- *                                     // actually execute the refinement,
- * tria->execute_coarsening_and_refinement ();
- *                                     // redistribute dofs,
- * dof_handler->distribute_dofs (fe);
- *                                     // and interpolate the solution
- * Vector<double> interpolated_solution(dof_handler->n_dofs());
- * soltrans.interpolate(interpolated_solution);
- * @endverbatim
+ SolutionTransfer<dim, VECTOR> soltrans(dof_handler);
+                                     // flag some cells for refinement
+                                     // and coarsening, e.g.
+ GridRefinement::refine_and_coarsen_fixed_fraction(
+   tria, error_indicators, 0.3, 0.05);
+                                     // prepare the triangulation,
+ tria.prepare_coarsening_and_refinement();
+                                     // prepare the SolutionTransfer object
+                                     // for coarsening and refinement and give
+                                     // the solution vector that we intend to
+                                     // interpolate later,
+ soltrans.prepare_for_coarsening_and_refinement(solution);
+                                     // actually execute the refinement,
+ tria.execute_coarsening_and_refinement ();
+                                     // redistribute dofs,
+ dof_handler.distribute_dofs (fe);
+                                     // and interpolate the solution
+ VECTOR interpolated_solution;
+ //create VECTOR in the right size here
+ soltrans.interpolate(interpolated_solution);
+@endverbatim
+ *
+ * <h3>Use for Serialization</h3>
+ *
+ * This class can be used to serialize and later deserialize a distributed mesh with solution vectors
+ * to a file. If you use more than one DoFHandler and therefore more than one SolutionTransfer object, they need to be serialized and deserialized in the same order.
+ *
+ * If vector has the locally relevant DoFs, serialization works as follows:
+ *@verbatim
+parallel::distributed::SolutionTransfer<dim,VECTOR> sol_trans(dof_handler);
+sol_trans.prepare_serialization (vector);
+
+triangulation.save(filename);
+@endverbatim
+* For deserialization the vector needs to be a distributed vector (without ghost elements):
+* @verbatim
+//[create coarse mesh...]
+triangulation.load(filename);
+
+parallel::distributed::SolutionTransfer<dim,VECTOR> sol_trans(dof_handler);
+sol_trans.deserialize (distributed_vector);
+@endverbatim
+ * 
  * @ingroup distributed
- * @author Timo Heister, 2009
+ * @author Timo Heister, 2009-2011
  */
     template<int dim, typename VECTOR, class DH=DoFHandler<dim> >
     class SolutionTransfer
@@ -113,12 +152,49 @@ namespace parallel
 
 
                                         /**
-                                         * return the size in bytes that need
+                                         * Return the size in bytes that need
                                          * to be stored per cell.
                                          */
        unsigned int get_data_size() const;
 
 
+                                        /**
+                                        * Prepare the serialization of the
+                                        * given vector. The serialization is
+                                        * done by Triangulation::save(). The
+                                        * given vector needs all information
+                                        * on the locally active DoFs (it must
+                                        * be ghosted). See documentation of
+                                        * this class for more information.
+                                        */
+       void prepare_serialization(const VECTOR &in);
+
+
+                                        /**
+                                         * Same as the function above, only
+                                         * for a list of vectors.
+                                         */
+       void prepare_serialization(const std::vector<const VECTOR*> &all_in);
+
+
+                                        /**
+                                         * Execute the deserialization of the
+                                         * given vector. This needs to be
+                                         * done after calling
+                                         * Triangulation::load(). The given
+                                         * vector must be a fully distributed
+                                         * vector without ghost elements. See
+                                         * documentation of this class for
+                                         * more information.
+                                         */
+       void deserialize(VECTOR &in);
+
+
+                                        /**
+                                         * Same as the function above, only
+                                         * for a list of vectors.
+                                         */
+       void deserialize(std::vector<VECTOR*> &all_in);
 
       private:
                                         /**
@@ -173,6 +249,12 @@ namespace parallel
                             const void* data,
                             std::vector<VECTOR*> &all_out);
 
+
+       /**
+        *
+        */
+       void register_data_attach(const std::size_t size);
+
     };
 
 
index 1d98cd845a71981ba4fc95e0a445270689984dec..1208278c45e9b964a8d7f21fc502db6eb92716fb 100644 (file)
@@ -436,6 +436,21 @@ namespace parallel
                                          */
        unsigned int get_checksum () const;
 
+                                        /**
+                                         * Save the refinement information from the coarse mesh into the given
+                                         * file. This file needs to be reachable from all nodes in the computation
+                                         * on a shared network file system. See the SolutionTransfer class
+                                         * on how to store solution vectors into this file.
+                                         */
+       void save(const char* filename) const;
+
+                                         /**
+                                          * Load the refinement information saved with save() back in. The mesh
+                                          * must contain the same coarse mesh that was used in save(). You need
+                                          * to load with the same number of CPUs that you saved with.
+                                          */
+       void load(const char* filename);
+       
                                         /**
                                          * Used to inform in the callbacks of
                                          * register_data_attach() and
index d1253a3fdb73a9f0c09e2186d9d4e2cccc74a501..0bcc2658de7b450944e035173f12b783035ce806 100644 (file)
@@ -93,10 +93,20 @@ namespace parallel
     SolutionTransfer<dim, VECTOR, DH>::
     prepare_for_coarsening_and_refinement (const std::vector<const VECTOR*> &all_in)
     {
-      Assert(all_in.size() > 0, ExcMessage("Please transfer atleast one vector!"));
+      
       input_vectors = all_in;
-      SolutionTransfer<dim, VECTOR, DH> *ptr = this;
+         register_data_attach( get_data_size() * input_vectors.size() );
+    }
+
+
 
+       template<int dim, typename VECTOR, class DH>
+    void
+    SolutionTransfer<dim, VECTOR, DH>::
+       register_data_attach(const std::size_t size)
+       {
+         Assert(size > 0, ExcMessage("Please transfer at least one vector!"));
+         
 //TODO: casting away constness is bad
       parallel::distributed::Triangulation<dim> * tria
        = (dynamic_cast<parallel::distributed::Triangulation<dim>*>
@@ -105,16 +115,16 @@ namespace parallel
       Assert (tria != 0, ExcInternalError());
 
       offset
-       = tria->register_data_attach(static_cast<size_t>
-                                    (get_data_size() * input_vectors.size()),
+       = tria->register_data_attach(size,
                                     std_cxx1x::bind(&SolutionTransfer<dim, VECTOR, DH>::pack_callback,
-                                                    ptr,
+                                                    this,
                                                     std_cxx1x::_1,
                                                     std_cxx1x::_2,
                                                     std_cxx1x::_3));
-    }
-
+         
+       }
 
+  
 
     template<int dim, typename VECTOR, class DH>
     void
@@ -127,6 +137,52 @@ namespace parallel
 
 
 
+    template<int dim, typename VECTOR, class DH>
+    void
+    SolutionTransfer<dim, VECTOR, DH>::
+       prepare_serialization(const VECTOR &in)
+       {
+      std::vector<const VECTOR*> all_in(1, &in);
+      prepare_serialization(all_in);
+       }
+
+
+
+       template<int dim, typename VECTOR, class DH>
+    void
+    SolutionTransfer<dim, VECTOR, DH>::
+       prepare_serialization(const std::vector<const VECTOR*> &all_in)
+       {
+         prepare_for_coarsening_and_refinement (all_in);
+       }
+
+
+
+    template<int dim, typename VECTOR, class DH>
+    void
+    SolutionTransfer<dim, VECTOR, DH>::
+       deserialize(VECTOR &in)
+       {
+      std::vector<VECTOR*> all_in(1, &in);
+      deserialize(all_in);
+       }
+
+
+
+       template<int dim, typename VECTOR, class DH>
+    void
+    SolutionTransfer<dim, VECTOR, DH>::
+       deserialize(std::vector<VECTOR*> &all_in)
+       {
+         register_data_attach( get_data_size() * all_in.size() );
+
+         // this makes interpolate() happy
+         input_vectors.resize(all_in.size());
+         
+         interpolate(all_in);
+       }
+
+
     template<int dim, typename VECTOR, class DH>
     void
     SolutionTransfer<dim, VECTOR, DH>::
index de7b273f7f4d9d53b33893554e6291b4424f411b..63db8611f36b59def2a0920fd5ee5a9850dec445 100644 (file)
@@ -37,6 +37,7 @@
 #include <algorithm>
 #include <numeric>
 #include <iostream>
+#include <fstream>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -1947,7 +1948,92 @@ namespace parallel
     }
 
 
+    template <int dim, int spacedim>
+    void
+    Triangulation<dim,spacedim>::
+    save(const char* filename) const
+    {
+         int real_data_size = 0;
+         if (attached_data_size>0)
+               real_data_size = attached_data_size+sizeof(CellStatus);
+         
+         if (my_subdomain==0)
+         {
+               std::string fname=std::string(filename)+".info";
+               std::ofstream f(fname.c_str());
+               f << Utilities::System::get_n_mpi_processes (mpi_communicator) << " "
+               << real_data_size << " "
+               << attached_data_pack_callbacks.size() << std::endl;
+         }
+
+         if (attached_data_size>0)
+         {
+               const_cast<dealii::parallel::distributed::Triangulation<dim, spacedim>*>(this)
+               ->attach_mesh_data();
+         }
+
+      dealii::internal::p4est::functions<dim>::save(filename, parallel_forest, attached_data_size>0);
+
+         dealii::parallel::distributed::Triangulation<dim, spacedim>* tria
+               = const_cast<dealii::parallel::distributed::Triangulation<dim, spacedim>*>(this);
+         
+         tria->n_attached_datas = 0;
+         tria->attached_data_size = 0;
+         tria->attached_data_pack_callbacks.clear();
+         
+                                          // and release the data
+         void * userptr = parallel_forest->user_pointer;
+         dealii::internal::p4est::functions<dim>::reset_data (parallel_forest, 0, NULL, NULL);
+         parallel_forest->user_pointer = userptr;
+    }
+
+
+    template <int dim, int spacedim>
+    void
+    Triangulation<dim,spacedim>::
+    load(const char* filename)
+    {
+      dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
+      parallel_forest = 0;
+      dealii::internal::p4est::functions<dim>::connectivity_destroy (connectivity);
+      connectivity=0;
+      
+         unsigned int numcpus, attached_size, attached_count;
+         {
+               std::string fname=std::string(filename)+".info";
+               std::ifstream f(fname.c_str());
+               f >> numcpus >> attached_size >> attached_count;
+               if (numcpus != Utilities::System::get_n_mpi_processes (mpi_communicator))
+                 throw ExcInternalError();
+         }
+
+         attached_data_size = 0;
+         n_attached_datas = 0;
+         
+      parallel_forest = dealii::internal::p4est::functions<dim>::load (
+               filename, mpi_communicator,
+               attached_size, attached_size>0,
+               this,
+               &connectivity);
+
+         try
+         {
+               copy_local_forest_to_triangulation ();
+         }
+      catch (const typename Triangulation<dim>::DistortedCellList &)
+         {
+                                          // the underlying
+                                          // triangulation should not
+                                          // be checking for
+                                          // distorted cells
+               AssertThrow (false, ExcInternalError());
+         }
+
+      update_number_cache ();
+    }
+
 
+    
     template <int dim, int spacedim>
     unsigned int
     Triangulation<dim,spacedim>::get_checksum () const

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.