[cig-commits] commit 1880 by dannberg to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Sep 11 07:57:10 PDT 2013
Revision 1880
add the mesh refinement criterion nonadiabatic temperature and change the compositin mesh refinement criterion so that it can use several fields
U branches/j-dannberg/include/aspect/mesh_refinement/composition.h
U branches/j-dannberg/include/aspect/mesh_refinement/interface.h
A branches/j-dannberg/include/aspect/mesh_refinement/nonadiabatic_temperature.h
U branches/j-dannberg/source/mesh_refinement/composition.cc
U branches/j-dannberg/source/mesh_refinement/interface.cc
A branches/j-dannberg/source/mesh_refinement/nonadiabatic_temperature.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1880&peg=1880
Diff:
Modified: branches/j-dannberg/include/aspect/mesh_refinement/composition.h
===================================================================
--- branches/j-dannberg/include/aspect/mesh_refinement/composition.h 2013-09-11 14:48:13 UTC (rev 1879)
+++ branches/j-dannberg/include/aspect/mesh_refinement/composition.h 2013-09-11 14:55:19 UTC (rev 1880)
@@ -41,6 +41,13 @@
class Composition : public Interface<dim>,
public SimulatorAccess<dim>
{
+ /**
+ * Initialization function.
+ */
+ virtual
+ void
+ initialize ();
+
public:
/**
* Execute this mesh refinement criterion.
@@ -54,6 +61,28 @@
virtual
void
execute (Vector<float> &error_indicators) 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:
+ /**
+ * The scaling factors that should be applied to the individual
+ * refinement indicators before merging.
+ */
+ std::vector<double> composition_scaling_factors;
};
}
}
Modified: branches/j-dannberg/include/aspect/mesh_refinement/interface.h
===================================================================
--- branches/j-dannberg/include/aspect/mesh_refinement/interface.h 2013-09-11 14:48:13 UTC (rev 1879)
+++ branches/j-dannberg/include/aspect/mesh_refinement/interface.h 2013-09-11 14:55:19 UTC (rev 1880)
@@ -72,6 +72,13 @@
~Interface ();
/**
+ * Initialization function.
+ */
+ virtual
+ void
+ initialize ();
+
+ /**
* Execute this mesh refinement criterion.
*
* @param[out] error_indicators A vector that for every active
Added: branches/j-dannberg/include/aspect/mesh_refinement/nonadiabatic_temperature.h
===================================================================
--- branches/j-dannberg/include/aspect/mesh_refinement/nonadiabatic_temperature.h (rev 0)
+++ branches/j-dannberg/include/aspect/mesh_refinement/nonadiabatic_temperature.h 2013-09-11 14:55:19 UTC (rev 1880)
@@ -0,0 +1,62 @@
+/*
+ Copyright (C) 2013 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: temperature.h 1433 2012-12-08 08:24:55Z bangerth $ */
+
+
+#ifndef __aspect__mesh_refinement_nonadiabatic_temperature_h
+#define __aspect__mesh_refinement_nonadiabatic_temperature_h
+
+#include <aspect/mesh_refinement/interface.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+
+ /**
+ * A class that implements a mesh refinement criterion based on
+ * the excess temperature (difference between temperature and
+ * adiabatic temperature).
+ *
+ * @ingroup MeshRefinement
+ */
+ template <int dim>
+ class NonadiabaticTemperature : public Interface<dim>,
+ public SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * 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) const;
+ };
+ }
+}
+
+#endif
Modified: branches/j-dannberg/source/mesh_refinement/composition.cc
===================================================================
--- branches/j-dannberg/source/mesh_refinement/composition.cc 2013-09-11 14:48:13 UTC (rev 1879)
+++ branches/j-dannberg/source/mesh_refinement/composition.cc 2013-09-11 14:55:19 UTC (rev 1880)
@@ -31,6 +31,19 @@
{
template <int dim>
void
+ Composition<dim>::initialize ()
+ {
+ AssertThrow (composition_scaling_factors.size() == this->n_compositional_fields()
+ ||
+ composition_scaling_factors.size() == 0,
+ ExcMessage ("The number of scaling factors given here must either be "
+ "zero or equal to the number of chosen refinement criteria."));
+ if (composition_scaling_factors.size() == 0)
+ composition_scaling_factors = std::vector<double> (this->n_compositional_fields(), 1.0);
+ }
+
+ template <int dim>
+ void
Composition<dim>::execute(Vector<float> &indicators) const
{
AssertThrow (this->n_compositional_fields() >= 1,
@@ -52,9 +65,57 @@
0,
0,
this->get_triangulation().locally_owned_subdomain());
+ for (unsigned int i=0; i<indicators.size(); ++i)
+ this_indicator[i] *= composition_scaling_factors[c];
indicators += this_indicator;
}
}
+
+ template <int dim>
+ void
+ Composition<dim>::
+ declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Mesh refinement");
+ {
+ prm.enter_subsection("Composition");
+ {
+ prm.declare_entry("Compositional field scaling factors",
+ "",
+ Patterns::List (Patterns::Double(0)),
+ "A list of scaling factors by which every individual compositional "
+ "field will be multiplied by. If only a single compositional "
+ "field exists, then this parameter has no particular meaning. "
+ "On the other hand, if multiple criteria are chosen, then these "
+ "factors are used to weigh the various indicators relative to "
+ "each other. "
+ "
"
+ "If the list of indicators given in this parameter is empty, then this "
+ "indicates that they should all be chosen equal to one. If the list "
+ "is not empty then it needs to have as many entries as there are "
+ "compositional fields.");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+ template <int dim>
+ void
+ Composition<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Mesh refinement");
+ {
+ prm.enter_subsection("Composition");
+ {
+ composition_scaling_factors
+ = Utilities::string_to_double(
+ Utilities::split_string_list(prm.get("Compositional field scaling factors")));
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
}
}
Modified: branches/j-dannberg/source/mesh_refinement/interface.cc
===================================================================
--- branches/j-dannberg/source/mesh_refinement/interface.cc 2013-09-11 14:48:13 UTC (rev 1879)
+++ branches/j-dannberg/source/mesh_refinement/interface.cc 2013-09-11 14:55:19 UTC (rev 1880)
@@ -22,6 +22,7 @@
#include <aspect/mesh_refinement/interface.h>
#include <aspect/simulator_access.h>
+#include <aspect/mesh_refinement/composition.h>
#include <typeinfo>
@@ -39,6 +40,12 @@
template <int dim>
void
+ Interface<dim>::initialize ()
+ {}
+
+
+ template <int dim>
+ void
Interface<dim>::declare_parameters (ParameterHandler &)
{}
@@ -65,8 +72,11 @@
for (typename std::list<std_cxx1x::shared_ptr<Interface<dim> > >::iterator
p = mesh_refinement_objects.begin();
p != mesh_refinement_objects.end(); ++p)
+ {
if (dynamic_cast<const SimulatorAccess<dim>*>(p->get()) != 0)
dynamic_cast<SimulatorAccess<dim>&>(**p).initialize (simulator);
+ (*p)->initialize();
+ }
// also extract the MPI communicator over which this simulator
// operates so that we can later scale vectors that live on
@@ -317,6 +327,7 @@
// find out which plugins are requested and the various other
// parameters we declare here
std::vector<std::string> plugin_names;
+ std::vector<int> compositional_fields;
prm.enter_subsection("Mesh refinement");
{
plugin_names
@@ -349,11 +360,13 @@
AssertThrow (plugin_names.size() >= 1,
ExcMessage ("You need to provide at least one mesh refinement criterion in the input file!"));
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],
"Mesh refinement::Refinement criteria merge operation",
prm)));
+ }
}
Added: branches/j-dannberg/source/mesh_refinement/nonadiabatic_temperature.cc
===================================================================
--- branches/j-dannberg/source/mesh_refinement/nonadiabatic_temperature.cc (rev 0)
+++ branches/j-dannberg/source/mesh_refinement/nonadiabatic_temperature.cc 2013-09-11 14:55:19 UTC (rev 1880)
@@ -0,0 +1,137 @@
+/*
+ Copyright (C) 2013 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: temperature.cc 1437 2012-12-08 12:02:49Z bangerth $ */
+
+
+#include <aspect/mesh_refinement/nonadiabatic_temperature.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/derivative_approximation.h>
+
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+ template <int dim>
+ void
+ NonadiabaticTemperature<dim>::execute(Vector<float> &indicators) const
+ {
+ indicators = 0;
+
+ // create a vector in which we set the temperature block to
+ // be a finite element interpolation of the density.
+ // we do so by setting up a quadrature formula with the
+ // temperature unit support points, then looping over these
+ // points, compute the output quantity at them, and writing
+ // the result into the output vector in the same order
+ // (because quadrature points and temperature dofs are,
+ // by design of the quadrature formula, numbered in the
+ // same way)
+ LinearAlgebra::BlockVector vec_distributed (this->introspection().index_sets.system_partitioning,
+ this->get_mpi_communicator());
+
+ const Quadrature<dim> quadrature(this->get_fe().base_element(2).get_unit_support_points());
+ std::vector<unsigned int> local_dof_indices (this->get_fe().dofs_per_cell);
+ FEValues<dim> fe_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature,
+ update_quadrature_points | update_values);
+
+ typename MaterialModel::Interface<dim>::MaterialModelInputs in(quadrature.size(),
+ this->n_compositional_fields());
+
+ // the values of the compositional fields are stored as blockvectors for each field
+ // we have to extract them in this structure
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = this->get_dof_handler().begin_active(),
+ endc = this->get_dof_handler().end();
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit(cell);
+ fe_values[this->introspection().extractors.temperature].get_function_values (this->get_solution(),
+ in.temperature);
+ in.position = fe_values.get_quadrature_points();
+ in.strain_rate.resize(0);// we are not reading the viscosity
+
+ cell->get_dof_indices (local_dof_indices);
+
+ // for each temperature dof, write into the output
+ // vector the density. note that quadrature points and
+ // dofs are enumerated in the same order
+ for (unsigned int i=0; i<this->get_fe().base_element(2).dofs_per_cell; ++i)
+ {
+ const unsigned int system_local_dof
+ = this->get_fe().component_to_system_index(/*temperature component=*/dim+1,
+ /*dof index within component=*/i);
+
+ vec_distributed(local_dof_indices[system_local_dof])
+ = in.temperature[i] - this->get_adiabatic_conditions().temperature(in.position[i]);
+ }
+ }
+
+ // now create a vector with the requisite ghost elements
+ // and use it for estimating the gradients
+ LinearAlgebra::BlockVector vec (this->introspection().index_sets.system_relevant_partitioning,
+ this->get_mpi_communicator());
+ vec = vec_distributed;
+
+ DerivativeApproximation::approximate_gradient (this->get_mapping(),
+ this->get_dof_handler(),
+ vec,
+ indicators,
+ //TODO: replace by the appropriate component mask
+ dim+1);
+
+ // Scale gradient in each cell with the correct power of h. Otherwise,
+ // error indicators do not reduce when refined if there is a density
+ // jump. We need at least order 1 for the error not to grow when
+ // refining, so anything >1 should work. (note that the gradient
+ // itself scales like 1/h, so multiplying it with any factor h^s, s>1
+ // will yield convergence of the error indicators to zero as h->0)
+ const double power = 1.0 + dim/2.0;
+ {
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = this->get_dof_handler().begin_active(),
+ endc = this->get_dof_handler().end();
+ unsigned int i=0;
+ for (; cell!=endc; ++cell, ++i)
+ if (cell->is_locally_owned())
+ indicators(i) *= std::pow(cell->diameter(), power);
+ }
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+ ASPECT_REGISTER_MESH_REFINEMENT_CRITERION(NonadiabaticTemperature,
+ "nonadiabatic temperature",
+ "A mesh refinement criterion that computes "
+ "refinement indicators from the excess temperature"
+ "(difference between temperature and adiabatic "
+ "temperature.")
+ }
+}
More information about the CIG-COMMITS
mailing list