[cig-commits] [commit] master: add a minimum refinement function mesh refinement criterium (2ce4f2a)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue May 20 07:34:33 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/426fbfcad3568313f61490f77127c5e1c97a005d...a871ef437febf4fd67ce6f1ff7349b8a57987c8a
>---------------------------------------------------------------
commit 2ce4f2a57686030d2571f24b4f9f3df6c10f2eeb
Author: Juliane Dannberg <dannberg at gfz-potsdam.de>
Date: Mon May 19 21:22:35 2014 +0200
add a minimum refinement function mesh refinement criterium
>---------------------------------------------------------------
2ce4f2a57686030d2571f24b4f9f3df6c10f2eeb
include/aspect/mesh_refinement/interface.h | 22 +++-
.../minimum_refinement_function.h} | 54 +++++-----
source/mesh_refinement/interface.cc | 82 +++++++++++++++
.../mesh_refinement/minimum_refinement_function.cc | 112 +++++++++++++++++++++
source/simulator/core.cc | 3 +
5 files changed, 245 insertions(+), 28 deletions(-)
diff --git a/include/aspect/mesh_refinement/interface.h b/include/aspect/mesh_refinement/interface.h
index 895886d..81b5f46 100644
--- a/include/aspect/mesh_refinement/interface.h
+++ b/include/aspect/mesh_refinement/interface.h
@@ -82,7 +82,18 @@ namespace aspect
*/
virtual
void
- execute (Vector<float> &error_indicators) const = 0;
+ execute (Vector<float> &error_indicators) const;
+
+ /**
+ * After cells have been marked for coarsening/refinement, apply
+ * additional criteria independent of the error estimate.
+ *
+ * @param[in] max_grid_level The maximum refinement level of the
+ * mesh
+ */
+ virtual
+ void
+ tag_additional_cells (unsigned int max_grid_level) const;
/**
* Declare the parameters this class takes through input files.
@@ -154,6 +165,15 @@ namespace aspect
execute (Vector<float> &error_indicators) const;
/**
+ * Apply additional refinement criteria independent of the error
+ * estimate for all of the mesh refinement objects that have been
+ * requested in the input file.
+ */
+ virtual
+ void
+ tag_additional_cells (unsigned int max_grid_level) const;
+
+ /**
* Declare the parameters of all known mesh refinement plugins, as
* well as of ones this class has itself.
*/
diff --git a/include/aspect/heating_model/constant_heating.h b/include/aspect/mesh_refinement/minimum_refinement_function.h
similarity index 53%
copy from include/aspect/heating_model/constant_heating.h
copy to include/aspect/mesh_refinement/minimum_refinement_function.h
index 8540bc7..21f645f 100644
--- a/include/aspect/heating_model/constant_heating.h
+++ b/include/aspect/mesh_refinement/minimum_refinement_function.h
@@ -17,46 +17,46 @@
along with ASPECT; see the file doc/COPYING. If not see
<http://www.gnu.org/licenses/>.
*/
-/* $Id$ */
+/* $Id: composition.h 1880 2013-09-11 14:55:19Z dannberg $ */
-#ifndef __aspect__heating_model_constant_heating_h
-#define __aspect__heating_model_constant_heating_h
+#ifndef __aspect__mesh_refinement_minimum_refinement_function_h
+#define __aspect__mesh_refinement_minimum_refinement_function_h
-#include <aspect/heating_model/interface.h>
+#include <aspect/mesh_refinement/interface.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/base/parsed_function.h>
namespace aspect
{
- namespace HeatingModel
+ namespace MeshRefinement
{
- using namespace dealii;
/**
- * A class that implements a constant radiogenic heating rate.
+ * A class that implements a minimum refinement level
+ * based on a functional description provided in the
+ * input file.
*
- * @ingroup HeatingModels
+ * @ingroup MeshRefinement
*/
template <int dim>
- class ConstantHeating : public Interface<dim>
+ class MinimumRefinementFunction : public Interface<dim>,
+ public SimulatorAccess<dim>
{
public:
/**
- * Return the specific heating rate. For the
- * current class, this function obviously simply returns a
- * constant value.
+ * After cells have been marked for coarsening/refinement, apply
+ * additional criteria independent of the error estimate.
+ *
+ * @param[in] max_grid_level The maximum refinement level of the
+ * mesh
*/
virtual
- double
- specific_heating_rate (const double,
- const double,
- const std::vector<double> &,
- const Point<dim> &) const;
+ void
+ tag_additional_cells (unsigned int max_grid_level) const;
/**
- * @name Functions used in dealing with run-time parameters
- * @{
- */
- /**
* Declare the parameters this class takes through input files.
*/
static
@@ -64,20 +64,20 @@ namespace aspect
declare_parameters (ParameterHandler &prm);
/**
- * Read the parameters this class declares from the parameter file.
+ * Read the parameters this class declares from the parameter
+ * file.
*/
virtual
void
parse_parameters (ParameterHandler &prm);
- /**
- * @}
- */
private:
- double radiogenic_heating_rate;
+ /**
+ * A function object representing the minimum refinement level.
+ */
+ Functions::ParsedFunction<1> min_refinement_level;
};
}
}
-
#endif
diff --git a/source/mesh_refinement/interface.cc b/source/mesh_refinement/interface.cc
index 4b2f70b..f2d653f 100644
--- a/source/mesh_refinement/interface.cc
+++ b/source/mesh_refinement/interface.cc
@@ -39,6 +39,21 @@ namespace aspect
template <int dim>
void
+ Interface<dim>::execute (Vector<float> &error_indicators) const
+ {
+ for (unsigned int i=0; i<error_indicators.size(); ++i)
+ error_indicators[i] = 0;
+ }
+
+
+ template <int dim>
+ void
+ Interface<dim>::tag_additional_cells (unsigned int) const
+ {}
+
+
+ template <int dim>
+ void
Interface<dim>::declare_parameters (ParameterHandler &)
{}
@@ -196,6 +211,73 @@ namespace aspect
}
+ template <int dim>
+ void
+ Manager<dim>::tag_additional_cells (unsigned int max_grid_level) const
+ {
+ Assert (mesh_refinement_objects.size() > 0, ExcInternalError());
+
+ // call the tag_additional_cells() functions of all
+ // plugins we have here in turns.
+ unsigned int index = 0;
+ for (typename std::list<std_cxx1x::shared_ptr<Interface<dim> > >::const_iterator
+ p = mesh_refinement_objects.begin();
+ p != mesh_refinement_objects.end(); ++p, ++index)
+ {
+ try
+ {
+ (*p)->tag_additional_cells (max_grid_level);
+ }
+
+ // 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);
+ }
+ }
+ }
+
+
+
+
+
// -------------------------------- Deal with registering plugins and automating
// -------------------------------- their setup and selection at run time
diff --git a/source/mesh_refinement/minimum_refinement_function.cc b/source/mesh_refinement/minimum_refinement_function.cc
new file mode 100644
index 0000000..0ba2180
--- /dev/null
+++ b/source/mesh_refinement/minimum_refinement_function.cc
@@ -0,0 +1,112 @@
+/*
+ Copyright (C) 2014 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: composition.cc 1880 2013-09-11 14:55:19Z dannberg $ */
+
+
+#include <aspect/mesh_refinement/minimum_refinement_function.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_values.h>
+
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+ template <int dim>
+ void
+ MinimumRefinementFunction<dim>::tag_additional_cells (unsigned int max_grid_level) const
+ {
+ // evaluate a single point per cell
+ const QMidpoint<dim> quadrature_formula;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FEValues<dim> fe_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature_formula,
+ update_values |
+ update_quadrature_points );
+
+ // ensure minimum refinement level
+ const unsigned int max_level = std::min(max_grid_level+1,this->get_triangulation().n_levels());
+ for (unsigned int level=0; level < max_level;++level)
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = this->get_triangulation().begin_active(level);
+ cell != this->get_triangulation().end_active(level); ++cell)
+ {
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit(cell);
+ const double depth = this->get_geometry_model().depth(fe_values.quadrature_point(0));
+ const Point<1> point(depth);
+ if (level <= std::rint(min_refinement_level.value(point)))
+ cell->clear_coarsen_flag ();
+ if (level < std::rint(min_refinement_level.value(point)))
+ cell->set_refine_flag ();
+ }
+ }
+ }
+
+ template <int dim>
+ void
+ MinimumRefinementFunction<dim>::
+ declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Mesh refinement");
+ {
+ prm.enter_subsection("Minimum refinement function");
+ {
+ Functions::ParsedFunction<1>::declare_parameters (prm, 1);
+ /* "The minimum refinement level each cell should have, "
+ "and that can not be exceeded by coarsening. " */
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+ template <int dim>
+ void
+ MinimumRefinementFunction<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Mesh refinement");
+ {
+ prm.enter_subsection("Minimum refinement function");
+ {
+ min_refinement_level.parse_parameters (prm);
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+ ASPECT_REGISTER_MESH_REFINEMENT_CRITERION(MinimumRefinementFunction,
+ "minimum refinement function",
+ "A mesh refinement criterion that ensures a "
+ "minimum refinement level described by an "
+ "explicit formula with the depth as argument.")
+ }
+}
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index e94a84b..1322804 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -998,6 +998,9 @@ namespace aspect
parameters.refinement_fraction,
parameters.coarsening_fraction);
+ mesh_refinement_manager.tag_additional_cells (max_grid_level);
+
+
// limit maximum refinement level
if (triangulation.n_levels() > max_grid_level)
for (typename Triangulation<dim>::active_cell_iterator
More information about the CIG-COMMITS
mailing list