[cig-commits] r1387 - in trunk/aspect: . include/aspect include/aspect/mesh_refinement source source/mesh_refinement
bangerth at dealii.org
bangerth at dealii.org
Wed Nov 28 06:04:46 PST 2012
Author: bangerth
Date: 2012-11-28 07:04:46 -0700 (Wed, 28 Nov 2012)
New Revision: 1387
Added:
trunk/aspect/include/aspect/mesh_refinement/
trunk/aspect/include/aspect/mesh_refinement/interface.h
trunk/aspect/source/mesh_refinement/
trunk/aspect/source/mesh_refinement/interface.cc
Modified:
trunk/aspect/Makefile
Log:
Start work on breaking mesh refinement up into a collection of plugins.
Modified: trunk/aspect/Makefile
===================================================================
--- trunk/aspect/Makefile 2012-11-28 07:35:08 UTC (rev 1386)
+++ trunk/aspect/Makefile 2012-11-28 14:04:46 UTC (rev 1387)
@@ -44,6 +44,7 @@
# list the directories and the various kinds of files
all-dirs := source \
source/simulator \
+ source/mesh_refinement \
source/geometry_model \
source/gravity_model \
source/boundary_temperature \
Copied: trunk/aspect/include/aspect/mesh_refinement/interface.h (from rev 1375, trunk/aspect/include/aspect/postprocess/interface.h)
===================================================================
--- trunk/aspect/include/aspect/mesh_refinement/interface.h (rev 0)
+++ trunk/aspect/include/aspect/mesh_refinement/interface.h 2012-11-28 14:04:46 UTC (rev 1387)
@@ -0,0 +1,249 @@
+/*
+ Copyright (C) 2011, 2012 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+/* $Id$ */
+
+
+#ifndef __aspect__mesh_refinement_interface_h
+#define __aspect__mesh_refinement_interface_h
+
+#include <aspect/global.h>
+#include <aspect/plugins.h>
+
+#include <deal.II/base/std_cxx1x/shared_ptr.h>
+#include <deal.II/base/table_handler.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/distributed/tria.h>
+
+
+namespace aspect
+{
+ using namespace dealii;
+
+ template <int dim> class Simulator;
+
+
+ /**
+ * A namespace for everything to do with the decision on how to refine
+ * the mesh every few time steps.
+ *
+ * @ingroup MeshRefinement
+ **/
+ namespace MeshRefinement
+ {
+
+ /**
+ * This class declares the public interface of mesh refinement plugins. These plugins
+ * must implement a function that can be called between time steps to refine
+ * the mesh based on the solution and/or the location of a cell.
+ *
+ * Access to the data of the simulator is granted by the @p protected member functions
+ * of the SimulatorAccessor class, i.e., classes implementing this interface will
+ * in general want to derive from both this Interface class as well as from the
+ * SimulatorAccess class.
+ *
+ * @ingroup MeshRefinement
+ */
+ template <int dim>
+ class Interface
+ {
+ public:
+ /**
+ * Destructor. Does nothing but is virtual so that derived classes
+ * destructors are also virtual.
+ **/
+ virtual
+ ~Interface ();
+
+ /**
+ * Execute this mesh refinement criterion.
+ *
+ * @param[out] error_indicators A vector that for every active
+ * cell of the current mesh
+ * (which may be a partition of a distributed mesh) provides an error
+ * indicator. This vector will already have the correct size when the
+ * function is called.
+ */
+ virtual
+ void
+ execute (Vector<float> &error_indicators) = 0;
+
+ /**
+ * Declare the parameters this class takes through input files.
+ * Derived classes should overload this function if they actually
+ * do take parameters; this class declares a fall-back function
+ * that does nothing, so that postprocessor classes that do not
+ * take any parameters do not have to do anything at all.
+ *
+ * This function is static (and needs to be static in derived
+ * classes) so that it can be called without creating actual
+ * objects (because declaring parameters happens before we read
+ * the input file and thus at a time when we don't even know yet
+ * which postprocessor objects we need).
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter
+ * file. The default implementation in this class does nothing,
+ * so that derived classes that do not need any parameters do
+ * not need to implement it.
+ */
+ virtual
+ void
+ parse_parameters (ParameterHandler &prm);
+ };
+
+
+
+
+
+
+ /**
+ * A class that manages all objects that provide functionality to refine meshes.
+ *
+ * @ingroup MeshRefinement
+ */
+ template <int dim>
+ class Manager
+ {
+ public:
+ /**
+ * Initialize the plugins handled by this object for a given simulator.
+ *
+ * @param simulator A reference to the main simulator object to which the
+ * postprocessor implemented in the derived class should be applied.
+ **/
+ void initialize (const Simulator<dim> &simulator);
+
+ /**
+ * Execute all of the mesh refinement objects that have been
+ * requested in the input file. The error indicators are then
+ * each individually normalized and merged according to the operation
+ * specified in the input file (e.g., via a plus, a maximum operation,
+ * etc).
+ */
+ virtual
+ void
+ execute (Vector<float> &error_indicators) = 0;
+
+ /**
+ * Declare the parameters of all known mesh refinement plugins, as
+ * well as of ones this class has itself.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter
+ * file. This determines which mesh refinement objects will be
+ * created; then let these objects read their parameters as
+ * well.
+ */
+ void
+ parse_parameters (ParameterHandler &prm);
+
+ /**
+ * A function that is used to register mesh refinement objects
+ * in such a way that the Manager can deal with all of them
+ * without having to know them by name. This allows the files
+ * in which individual plugins are implement to register
+ * these plugins, rather than also having to modify the
+ * Manager class by adding the new mesh refinement class.
+ *
+ * @param name The name under which this plugin is to
+ * be called in parameter files.
+ * @param description A text description of what this model
+ * does and that will be listed in the documentation of
+ * the parameter file.
+ * @param declare_parameters_function A pointer to a function
+ * that declares the parameters for this plugin.
+ * @param factory_function A pointer to a function that creates
+ * such a mesh refinement object and returns a pointer to it.
+ **/
+ static
+ void
+ register_mesh_refinement_criterion (const std::string &name,
+ const std::string &description,
+ void (*declare_parameters_function) (ParameterHandler &),
+ Interface<dim> *(*factory_function) ());
+
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcMeshRefinementNameNotFound,
+ std::string,
+ << "Could not find entry <"
+ << arg1
+ << "> among the names of registered mesh refinement objects.");
+ private:
+ /**
+ * An enum that describes the different ways in which we can
+ * merge the results of multiple mesh refinement criteria.
+ */
+ enum MergeOperation
+ { plus, max };
+
+ /**
+ * How to merge the results of multiple mesh refinement criteria.
+ */
+ MergeOperation merge_operation;
+
+ /**
+ * A list of mesh refinement objects that have been requested
+ * in the parameter file.
+ */
+ std::list<std_cxx1x::shared_ptr<Interface<dim> > > mesh_refinement_objects;
+
+ /**
+ * An MPI communicator that spans the set of processors on
+ * which the simulator object lives.
+ */
+ MPI_Comm mpi_communicator;
+ };
+
+
+
+ /**
+ * Given a class name, a name, and a description for the parameter file for a mesh
+ * refinement object, register it with
+ * the aspect::MeshRefinement::Manager class.
+ *
+ * @ingroup MeshRefinement
+ */
+#define ASPECT_REGISTER_MESH_REFINEMENT_CRITERION(classname,name,description) \
+ template class classname<2>; \
+ template class classname<3>; \
+ namespace ASPECT_REGISTER_MESH_REFINEMENT_CRITERION_ ## classname \
+ { \
+ aspect::internal::Plugins::RegisterHelper<Interface<2>,classname<2> > \
+ dummy_ ## classname ## _2d (&aspect::MeshRefinement::Manager<2>::register_mesh_refinement_criterion, \
+ name, description); \
+ aspect::internal::Plugins::RegisterHelper<Interface<3>,classname<3> > \
+ dummy_ ## classname ## _3d (&aspect::MeshRefinement::Manager<3>::register_mesh_refinement_criterion, \
+ name, description); \
+ }
+ }
+}
+
+
+#endif
Copied: trunk/aspect/source/mesh_refinement/interface.cc (from rev 1375, trunk/aspect/source/postprocess/interface.cc)
===================================================================
--- trunk/aspect/source/mesh_refinement/interface.cc (rev 0)
+++ trunk/aspect/source/mesh_refinement/interface.cc 2012-11-28 14:04:46 UTC (rev 1387)
@@ -0,0 +1,297 @@
+/*
+ Copyright (C) 2012 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+/* $Id$ */
+
+
+#include <aspect/mesh_refinement/interface.h>
+#include <aspect/simulator.h>
+
+#include <typeinfo>
+
+
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+// ------------------------------ Interface -----------------------------
+
+ template <int dim>
+ Interface<dim>::~Interface ()
+ {}
+
+
+ template <int dim>
+ void
+ Interface<dim>::declare_parameters (ParameterHandler &)
+ {}
+
+
+
+ template <int dim>
+ void
+ Interface<dim>::parse_parameters (ParameterHandler &)
+ {}
+
+
+
+// ------------------------------ Manager -----------------------------
+
+ template <int dim>
+ void
+ Manager<dim>::initialize (const Simulator<dim> &simulator)
+ {
+ for (typename std::list<std_cxx1x::shared_ptr<Interface<dim> > >::iterator
+ p = mesh_refinement_objects.begin();
+ p != mesh_refinement_objects.end(); ++p)
+ dynamic_cast<SimulatorAccess<dim>&>(**p).initialize (simulator);
+
+ // also extract the MPI communicator over which this simulator
+ // operates so that we can later scale vectors that live on
+ // different processors. getting at this communicator is a bit
+ // indirect but requires no extra interfaces (we need the
+ // intermediary MySimulatorAccess class since SimulatorAccess's
+ // get_mpi_communicator function is only protected; through the
+ // 'using' declaration we can promote it to a public member).
+ class MySimulatorAccess : public SimulatorAccess<dim>
+ {
+ public:
+ using SimulatorAccess<dim>::get_mpi_communicator;
+ };
+ MySimulatorAccess simulator_access;
+ simulator_access.initialize (simulator);
+ mpi_communicator = simulator_access.get_mpi_communicator();
+ }
+
+
+
+ template <int dim>
+ void
+ Manager<dim>::execute (Vector<float> &error_indicators)
+ {
+ // call the execute() functions of all plugins we have
+ // here in turns. then normalize the output vector and
+ // verify that its values are non-negative numbers
+ std::vector<Vector<float> > all_error_indicators (mesh_refinement_objects.size(),
+ Vector<float>(error_indicators.size()));
+ unsigned int index = 0;
+ for (typename std::list<std_cxx1x::shared_ptr<Interface<dim> > >::iterator
+ p = mesh_refinement_objects.begin();
+ p != mesh_refinement_objects.end(); ++p, ++index)
+ {
+ try
+ {
+ (*p)->execute (all_error_indicators[index]);
+
+ for (unsigned int i=0; i<error_indicators.size(); ++i)
+ Assert (all_error_indicators[index](i) >= 0,
+ ExcMessage ("Error indicators must be non-negative numbers!"));
+
+ const double global_max = Utilities::MPI::max (all_error_indicators[index].linfty_norm(),
+ mpi_communicator);
+ if (global_max != 0)
+ all_error_indicators[index] /= global_max;
+ }
+ // plugins that throw exceptions usually do not result in
+ // anything good because they result in an unwinding of the stack
+ // and, if only one processor triggers an exception, the
+ // destruction of objects often causes a deadlock. thus, if
+ // an exception is generated, catch it, print an error message,
+ // and abort the program
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on MPI process <"
+ << Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)
+ << "> while running mesh refinement plugin <"
+ << typeid(**p).name()
+ << ">: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ // terminate the program!
+ MPI_Abort (MPI_COMM_WORLD, 1);
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on MPI process <"
+ << Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)
+ << "> while running mesh refinement plugin <"
+ << typeid(**p).name()
+ << ">: " << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ // terminate the program!
+ MPI_Abort (MPI_COMM_WORLD, 1);
+ }
+ }
+
+ // now merge the results
+ }
+
+
+// -------------------------------- Deal with registering plugins and automating
+// -------------------------------- their setup and selection at run time
+
+ namespace
+ {
+ std_cxx1x::tuple
+ <void *,
+ void *,
+ internal::Plugins::PluginList<Interface<2> >,
+ internal::Plugins::PluginList<Interface<3> > > registered_plugins;
+ }
+
+
+
+ template <int dim>
+ void
+ Manager<dim>::declare_parameters (ParameterHandler &prm)
+ {
+ // first declare the postprocessors we know about to
+ // choose from
+ prm.enter_subsection("Mesh refinement");
+ {
+ // construct a string for Patterns::MultipleSelection that
+ // contains the names of all registered plugins
+ const std::string pattern_of_names
+ = std_cxx1x::get<dim>(registered_plugins).get_pattern_of_names (true);
+ prm.declare_entry("List of mesh refinement criteria",
+ "plus",
+ Patterns::MultipleSelection(pattern_of_names),
+ "A comma separated list of mesh refinement criteria that "
+ "will be run whenever mesh refinement is required. The "
+ "results of each of these criteria will, i.e., the refinement "
+ "indicators they produce for all the cells of the mesh "
+ "will then be normalized to a range between zero and one "
+ "and the results of different criteria will then be "
+ "merged through the operation selected in this section.\n\n"
+ "The following criteria are available:\n\n"
+ +
+ std_cxx1x::get<dim>(registered_plugins).get_description_string());
+
+ prm.declare_entry("Merge operation",
+ "plus",
+ Patterns::Selection("plus|max"),
+ "If multiple mesh refinement criteria are computed for each cell, "
+ "then one will have to decide which one should win when deciding "
+ "which cell to refine. The operation that selects from these competing "
+ "criteria is the one that is selected here. The options are:\n\n"
+ "\\begin{itemize}\n"
+ "\\item \\texttt{plus}: Add the various error indicators together and "
+ "refine those cells on which the sum of indicators is largest.\n"
+ "\\item \\texttt{max}: Take the maximum of the various error indicators and "
+ "refine those cells on which the maximal indicators is largest.\n"
+ "\\end{itemize}");
+ }
+ prm.leave_subsection();
+
+ // now declare the parameters of each of the registered
+ // plugins in turn
+ std_cxx1x::get<dim>(registered_plugins).declare_parameters (prm);
+ }
+
+
+
+ template <int dim>
+ void
+ Manager<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ Assert (std_cxx1x::get<dim>(registered_plugins).plugins != 0,
+ ExcMessage ("No mesh refinement plugins registered!?"));
+
+ // first find out which plugins are requested
+ std::vector<std::string> plugin_names;
+ prm.enter_subsection("Mesh refinement");
+ {
+ plugin_names
+ = Utilities::split_string_list(prm.get("List of mesh refinement criteria"));
+
+ // while we're in this section, also read the merge operation
+ if (prm.get("Merge operation") == "plus")
+ merge_operation = plus;
+ else if (prm.get("Merge operation") == "max")
+ merge_operation = max;
+ else
+ Assert (false, ExcNotImplemented());
+ }
+ prm.leave_subsection();
+
+ // go through the list, create objects and let them parse
+ // their own parameters
+ for (unsigned int name=0; name<plugin_names.size(); ++name)
+ mesh_refinement_objects.push_back (std_cxx1x::shared_ptr<Interface<dim> >
+ (std_cxx1x::get<dim>(registered_plugins)
+ .create_plugin (plugin_names[name],
+ prm)));
+ }
+
+
+ template <int dim>
+ void
+ Manager<dim>::register_mesh_refinement_criterion (const std::string &name,
+ const std::string &description,
+ void (*declare_parameters_function) (ParameterHandler &),
+ Interface<dim> *(*factory_function) ())
+ {
+ std_cxx1x::get<dim>(registered_plugins).register_plugin (name,
+ description,
+ declare_parameters_function,
+ factory_function);
+ }
+
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace internal
+ {
+ namespace Plugins
+ {
+ template <>
+ std::list<internal::Plugins::PluginList<MeshRefinement::Interface<2> >::PluginInfo> *
+ internal::Plugins::PluginList<MeshRefinement::Interface<2> >::plugins = 0;
+ template <>
+ std::list<internal::Plugins::PluginList<MeshRefinement::Interface<3> >::PluginInfo> *
+ internal::Plugins::PluginList<MeshRefinement::Interface<3> >::plugins = 0;
+ }
+ }
+
+ namespace MeshRefinement
+ {
+#define INSTANTIATE(dim) \
+ template class Interface<dim>; \
+ template class Manager<dim>;
+
+ ASPECT_INSTANTIATE(INSTANTIATE)
+ }
+}
More information about the CIG-COMMITS
mailing list