[cig-commits] commit 2569 by dannberg to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Thu May 8 11:15:51 PDT 2014
Revision 2569
add a mesh refinement criterion that ensures a minimum refinement level in dependence of depth, plus some bug fixes
D branches/j-dannberg/Makefile
U branches/j-dannberg/include/aspect/mesh_refinement/interface.h
A branches/j-dannberg/include/aspect/mesh_refinement/minimum_refinement_function.h
U branches/j-dannberg/source/mesh_refinement/interface.cc
A branches/j-dannberg/source/mesh_refinement/minimum_refinement_function.cc
U branches/j-dannberg/source/simulator/core.cc
U branches/j-dannberg/source/simulator/initial_conditions.cc
D branches/j-dannberg/tests/Makefile
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2569&peg=2569
Diff:
Modified: branches/j-dannberg/include/aspect/mesh_refinement/interface.h
===================================================================
--- branches/j-dannberg/include/aspect/mesh_refinement/interface.h 2014-05-06 04:07:59 UTC (rev 2568)
+++ branches/j-dannberg/include/aspect/mesh_refinement/interface.h 2014-05-08 18:15:48 UTC (rev 2569)
@@ -89,9 +89,20 @@
*/
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.
* Derived classes should overload this function if they actually
* do take parameters; this class declares a fall-back function
@@ -158,6 +169,15 @@
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.
*/
Added: branches/j-dannberg/include/aspect/mesh_refinement/minimum_refinement_function.h
===================================================================
--- branches/j-dannberg/include/aspect/mesh_refinement/minimum_refinement_function.h (rev 0)
+++ branches/j-dannberg/include/aspect/mesh_refinement/minimum_refinement_function.h 2014-05-08 18:15:48 UTC (rev 2569)
@@ -0,0 +1,83 @@
+/*
+ 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.h 1880 2013-09-11 14:55:19Z dannberg $ */
+
+
+#ifndef __aspect__mesh_refinement_minimum_refinement_function_h
+#define __aspect__mesh_refinement_minimum_refinement_function_h
+
+#include <aspect/mesh_refinement/interface.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/base/parsed_function.h>
+
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+
+ /**
+ * A class that implements a minimum refinement level
+ * based on a functional description provided in the
+ * input file.
+ *
+ * @ingroup MeshRefinement
+ */
+ template <int dim>
+ class MinimumRefinementFunction : public Interface<dim>,
+ public SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * 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.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter
+ * file.
+ */
+ virtual
+ void
+ parse_parameters (ParameterHandler &prm);
+
+ private:
+ /**
+ * A function object representing the minimum refinement level.
+ */
+ Functions::ParsedFunction<1> min_refinement_level;
+ };
+ }
+}
+
+#endif
Modified: branches/j-dannberg/source/mesh_refinement/interface.cc
===================================================================
--- branches/j-dannberg/source/mesh_refinement/interface.cc 2014-05-06 04:07:59 UTC (rev 2568)
+++ branches/j-dannberg/source/mesh_refinement/interface.cc 2014-05-08 18:15:48 UTC (rev 2569)
@@ -46,6 +46,21 @@
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 &)
{}
@@ -206,6 +221,73 @@
}
+ 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
Added: branches/j-dannberg/source/mesh_refinement/minimum_refinement_function.cc
===================================================================
--- branches/j-dannberg/source/mesh_refinement/minimum_refinement_function.cc (rev 0)
+++ branches/j-dannberg/source/mesh_refinement/minimum_refinement_function.cc 2014-05-08 18:15:48 UTC (rev 2569)
@@ -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.")
+ }
+}
Modified: branches/j-dannberg/source/simulator/core.cc
===================================================================
--- branches/j-dannberg/source/simulator/core.cc 2014-05-06 04:07:59 UTC (rev 2568)
+++ branches/j-dannberg/source/simulator/core.cc 2014-05-08 18:15:48 UTC (rev 2569)
@@ -484,6 +484,9 @@
if (parameters.n_compositional_fields > 0)
statistics.add_value("Number of composition degrees of freedom",
introspection.system_dofs_per_block[3]);
+ if (parameters.include_melt_transport)
+ statistics.add_value("Number of porosity degrees of freedom",
+ introspection.system_dofs_per_block[(parameters.n_compositional_fields > 0 ? 4 : 3)]);
// then interpolate the current boundary velocities. this adds to
@@ -501,7 +504,7 @@
{
p->second->set_current_time (time);
VectorFunctionFromVelocityFunctionObject<dim> vel
- (parameters.n_compositional_fields,
+ (parameters.n_compositional_fields+(parameters.include_melt_transport ? 1 : 0),
std_cxx1x::bind (&VelocityBoundaryConditions::Interface<dim>::boundary_velocity,
p->second,
std_cxx1x::_1));
@@ -850,8 +853,9 @@
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
introspection.component_masks.compositional_fields
.push_back (finite_element.component_mask(introspection.extractors.compositional_fields[c]));
- introspection.component_masks.porosity
- = finite_element.component_mask (introspection.extractors.porosity);
+ if (parameters.include_melt_transport)
+ introspection.component_masks.porosity
+ = finite_element.component_mask (introspection.extractors.porosity);
}
@@ -868,7 +872,9 @@
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
n_C[c] = introspection.system_dofs_per_block
[introspection.block_indices.compositional_fields[c]];
- const types::global_dof_index n_Phi = introspection.system_dofs_per_block[3+parameters.n_compositional_fields];
+ types::global_dof_index n_Phi;
+ if (parameters.include_melt_transport)
+ n_Phi = introspection.system_dofs_per_block[3+parameters.n_compositional_fields];
IndexSet system_index_set = dof_handler.locally_owned_dofs();
@@ -899,11 +905,14 @@
n_u+n_p+n_T+n_C_so_far+n_C[c]));
n_C_so_far += n_C[c];
}
- introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u+n_p+n_T+n_C_so_far,
- n_u+n_p+n_T+n_C_so_far+n_Phi));
- introspection.index_sets.system_relevant_partitioning
- .push_back(introspection.index_sets.system_relevant_set.get_view(n_u+n_p+n_T+n_C_so_far,
- n_u+n_p+n_T+n_C_so_far+n_Phi));
+ if (parameters.include_melt_transport)
+ {
+ introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u+n_p+n_T+n_C_so_far,
+ n_u+n_p+n_T+n_C_so_far+n_Phi));
+ introspection.index_sets.system_relevant_partitioning
+ .push_back(introspection.index_sets.system_relevant_set.get_view(n_u+n_p+n_T+n_C_so_far,
+ n_u+n_p+n_T+n_C_so_far+n_Phi));
+ }
}
}
@@ -971,6 +980,9 @@
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
@@ -1064,13 +1076,16 @@
= solution.block(introspection.block_indices.compositional_fields[c]);
}
- assemble_advection_system (AdvectionField::porosity());
- build_advection_preconditioner(AdvectionField::porosity(),
- Phi_preconditioner);
- solve_advection(AdvectionField::porosity());
+ if (parameters.include_melt_transport)
+ {
+ assemble_advection_system (AdvectionField::porosity());
+ build_advection_preconditioner(AdvectionField::porosity(),
+ Phi_preconditioner);
+ solve_advection(AdvectionField::porosity());
- current_linearization_point.block(introspection.block_indices.porosity)
- = solution.block(introspection.block_indices.porosity);
+ current_linearization_point.block(introspection.block_indices.porosity)
+ = solution.block(introspection.block_indices.porosity);
+ }
// the Stokes matrix depends on the viscosity. if the viscosity
// depends on other solution variables, then after we need to
Modified: branches/j-dannberg/source/simulator/initial_conditions.cc
===================================================================
--- branches/j-dannberg/source/simulator/initial_conditions.cc 2014-05-06 04:07:59 UTC (rev 2568)
+++ branches/j-dannberg/source/simulator/initial_conditions.cc 2014-05-08 18:15:48 UTC (rev 2569)
@@ -70,12 +70,14 @@
{
initial_solution.reinit(system_rhs,false);
- // base element in the finite element is 2 for temperature (n=0) and 3 for
- // compositional fields (n>0)
+ // base element in the finite element is 2 for temperature (n=0), 3 for
+ // compositional fields (n>0) and 4 for porosity (or 3, if there are no
+ // compositional fields)
//TODO: can we use introspection here, instead of the hard coded numbers?
- const unsigned int base_element = (n==0 ? 2 : (n==1+parameters.n_compositional_fields) ? 4 : 3);
+ const unsigned int base_element = (n==0 ? 2 : (n==1+parameters.n_compositional_fields
+ && parameters.n_compositional_fields>0) ? 4 : 3);
- // get the temperature/composition support points
+ // get the temperature/composition/porosity support points
const std::vector<Point<dim> > support_points
= finite_element.base_element(base_element).get_unit_support_points();
Assert (support_points.size() != 0,
@@ -94,8 +96,8 @@
{
fe_values.reinit (cell);
- // go through the temperature/composition dofs and set their global values
- // to the temperature/composition field interpolated at these points
+ // go through the temperature/composition/porosity dofs and set their global values
+ // to the temperature/composition/porosity field interpolated at these points
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<finite_element.base_element(base_element).dofs_per_cell; ++i)
{
@@ -106,9 +108,9 @@
double value =
(base_element == 2 ?
initial_conditions->initial_temperature(fe_values.quadrature_point(i))
- : (base_element == 4 ?
- porosity_initial_conditions->initial_porosity(fe_values.quadrature_point(i))
- : compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i),n-1)));
+ : (base_element == 3 && parameters.n_compositional_fields > 0 ?
+ compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i),n-1)
+ : porosity_initial_conditions->initial_porosity(fe_values.quadrature_point(i))));
initial_solution(local_dof_indices[system_local_dof]) = value;
if (base_element != 2)
More information about the CIG-COMMITS
mailing list