From 2f4273f81aaa816a011fd49eafbe89e118a90700 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Mon, 21 Oct 2024 15:55:08 +0200 Subject: [PATCH 1/7] Derive TractionManager from ManagerBase --- include/aspect/boundary_traction/interface.h | 81 ++++++++++---- include/aspect/plugins.h | 3 +- source/boundary_traction/ascii_data.cc | 13 ++- .../initial_lithostatic_pressure.cc | 12 +- source/boundary_traction/interface.cc | 104 +++++++++++------- source/simulator/assemblers/stokes.cc | 7 +- source/simulator/assembly.cc | 4 +- source/simulator/melt.cc | 9 +- source/simulator/newton.cc | 2 +- 9 files changed, 152 insertions(+), 83 deletions(-) diff --git a/include/aspect/boundary_traction/interface.h b/include/aspect/boundary_traction/interface.h index 68bff12a0be..85a8cdd13d6 100644 --- a/include/aspect/boundary_traction/interface.h +++ b/include/aspect/boundary_traction/interface.h @@ -71,32 +71,13 @@ namespace aspect }; template - class Manager : public SimulatorAccess + class Manager : public Plugins::ManagerBase>, public SimulatorAccess { public: - /** - * Destructor. Made virtual since this class has virtual member - * functions. - */ - ~Manager () override; - - /** - * A function that is called at the beginning of each time step and - * calls the corresponding functions of all created plugins. - * - * The point of this function is to allow complex boundary traction - * models to do an initialization step once at the beginning of each - * time step. An example would be a model that needs to call an - * external program to compute the traction change at a boundary. - */ - virtual - void - update (); - /** * A function that calls the boundary_traction functions of all the - * individual boundary traction objects and uses the stored operators - * to combine them. + * individual boundary traction objects that are active for boundary id + * @p boundary_indicator and uses the stored operators to combine them. */ Tensor<1,dim> boundary_traction (const types::boundary_id boundary_indicator, @@ -114,7 +95,11 @@ namespace aspect * If there are no prescribed boundary traction plugins * for a particular boundary, this boundary identifier will not appear * in the map. + * + * @deprecated: This function will be removed. Use the function + * get_active_plugin_names() of the base class ManagerBase instead. */ + DEAL_II_DEPRECATED const std::map>> & get_active_boundary_traction_names () const; @@ -126,10 +111,32 @@ namespace aspect * boundary models for this boundary. If there are no prescribed * boundary traction plugins for a particular boundary this boundary * identifier will not appear in the map. + * + * @deprecated: This function has been removed. Use the function + * get_active_plugins() of the base class ManagerBase instead. */ + DEAL_II_DEPRECATED const std::map>>> & get_active_boundary_traction_conditions () const; + /** + * Return a list of boundary indicators that indicate for + * each active plugin which boundary id + * it is responsible for. The list of active plugins can be + * requested by calling get_active_plugins(). + */ + const std::vector & + get_active_plugin_boundary_indicators() const; + + /** + * Return a list of component masks that indicate for + * each active plugin which components it is responsible for. + * The list of plugin objects can be + * requested by calling get_active_plugins(). + */ + const std::vector & + get_active_plugin_component_masks() const; + /** * Declare the parameters of all known boundary traction plugins, as * well as the ones this class has itself. @@ -144,7 +151,7 @@ namespace aspect * then let these objects read their parameters as well. */ void - parse_parameters (ParameterHandler &prm); + parse_parameters (ParameterHandler &prm) override; /** * For the current plugin subsystem, write a connection graph of all of the @@ -193,9 +200,32 @@ namespace aspect std::unique_ptr> (*factory_function) ()); private: + /** + * A list of boundary indicators that indicate for + * each plugin in the list of plugin_objects which boundary id + * it is responsible for. By default each plugin + * is active for all boundaries, but this list + * can be modified by derived classes to limit the application + * of plugins to specific boundaries. + */ + std::vector boundary_indicators; + + /** + * A list of boundary indicators that indicate for + * each plugin in the list of plugin_objects which components + * it is responsible for. By default each plugin + * is active for all components, but this list + * can be modified by derived classes to limit the application + * of plugins to specific boundaries. + */ + std::vector component_masks; + /** * A list of boundary traction objects that have been requested in the * parameter file. + * + * @deprecated: This variable is no longer used, but needed to issue a proper + * error message in the function get_active_boundary_traction_conditions(). */ std::map>>> boundary_traction_objects; @@ -206,6 +236,11 @@ namespace aspect * mapped to one of the plugins of traction boundary conditions (e.g. * "function"). If the components string is empty, it is assumed the * plugins are used for all components. + * + * @deprecated: Remove this variable when the deprecated functions + * get_active_boundary_traction_names and + * get_active_boundary_traction_conditions are removed. Use the base class + * variable plugin_names instead. */ std::map>> boundary_traction_indicators; diff --git a/include/aspect/plugins.h b/include/aspect/plugins.h index ec0e7fec888..da4d97125f4 100644 --- a/include/aspect/plugins.h +++ b/include/aspect/plugins.h @@ -26,11 +26,12 @@ #include #include -#include #include +#include #include +#include #include #include #include diff --git a/source/boundary_traction/ascii_data.cc b/source/boundary_traction/ascii_data.cc index 6d84b16ec7b..1cea4822749 100644 --- a/source/boundary_traction/ascii_data.cc +++ b/source/boundary_traction/ascii_data.cc @@ -35,13 +35,16 @@ namespace aspect void AsciiData::initialize () { - for (const auto &bv : this->get_boundary_traction_manager().get_active_boundary_traction_conditions()) + unsigned int i=0; + for (const auto &plugin : this->get_boundary_traction_manager().get_active_plugins()) { - for (const auto &plugin : bv.second) - if (plugin.get() == this) - boundary_ids.insert(bv.first); + if (plugin.get() == this) + boundary_ids.insert(this->get_boundary_traction_manager().get_active_plugin_boundary_indicators()[i]); + + ++i; } - AssertThrow(*(boundary_ids.begin()) != numbers::invalid_boundary_id, + + AssertThrow(boundary_ids.empty() == false, ExcMessage("Did not find the boundary indicator for the traction ascii data plugin.")); Utilities::AsciiDataBoundary::initialize(boundary_ids, diff --git a/source/boundary_traction/initial_lithostatic_pressure.cc b/source/boundary_traction/initial_lithostatic_pressure.cc index c03db50a8ac..5b5fc1c774b 100644 --- a/source/boundary_traction/initial_lithostatic_pressure.cc +++ b/source/boundary_traction/initial_lithostatic_pressure.cc @@ -49,13 +49,15 @@ namespace aspect // Ensure the initial lithostatic pressure traction boundary conditions are used, // and register for which boundary indicators these conditions are set. std::set traction_bi; - for (const auto &p : this->get_boundary_traction_manager().get_active_boundary_traction_conditions()) + unsigned int i=0; + for (const auto &plugin : this->get_boundary_traction_manager().get_active_plugins()) { - for (const auto &plugin : p.second) - if (plugin.get() == this) - traction_bi.insert(p.first); + if (plugin.get() == this) + traction_bi.insert(this->get_boundary_traction_manager().get_active_plugin_boundary_indicators()[i]); + + ++i; } - AssertThrow(*(traction_bi.begin()) != numbers::invalid_boundary_id, + AssertThrow(traction_bi.empty() == false, ExcMessage("Did not find any boundary indicators for the initial lithostatic pressure plugin.")); // Determine whether traction boundary conditions are only set on the bottom diff --git a/source/boundary_traction/interface.cc b/source/boundary_traction/interface.cc index 2e02b8da1e5..1917009295d 100644 --- a/source/boundary_traction/interface.cc +++ b/source/boundary_traction/interface.cc @@ -33,23 +33,6 @@ namespace aspect { namespace BoundaryTraction { - template - Manager::~Manager() - = default; - - - - template - void - Manager::update () - { - for (const auto &boundary : boundary_traction_objects) - for (const auto &p : boundary.second) - p->update(); - } - - - namespace { std::tuple @@ -78,20 +61,32 @@ namespace aspect const Point &position, const Tensor<1,dim> &normal_vector) const { - typename std::map>>>::const_iterator boundary_plugins = - boundary_traction_objects.find(boundary_indicator); + Tensor<1,dim> traction; + + bool found_plugin = false; + unsigned int i=0; + for (const auto &plugin: this->plugin_objects) + { + if (boundary_indicators[i] == boundary_indicator) + { + found_plugin = true; + const Tensor<1,dim> plugin_traction = plugin->boundary_traction(boundary_indicator, + position, + normal_vector); + for (unsigned int d=0; d traction = Tensor<1,dim>(); - - for (const auto &plugin : boundary_plugins->second) - traction += plugin->boundary_traction(boundary_indicator, - position,normal_vector); - return traction; } @@ -110,11 +105,32 @@ namespace aspect const std::map>>> & Manager::get_active_boundary_traction_conditions () const { + AssertThrow(false, ExcMessage("This function has been removed. Use the function " + "get_active_plugins() of the base class ManagerBase " + "instead.")); return boundary_traction_objects; } + template + const std::vector & + Manager::get_active_plugin_boundary_indicators() const + { + return boundary_indicators; + } + + + + template + const std::vector & + Manager::get_active_plugin_component_masks() const + { + return component_masks; + } + + + template void Manager::declare_parameters (ParameterHandler &prm) @@ -254,27 +270,39 @@ namespace aspect { boundary_traction_indicators[boundary_id] = std::make_pair(comp,std::vector(1,value)); } + + this->plugin_names.push_back(value); + boundary_indicators.push_back(boundary_id); + + const bool default_component_mask = comp.empty(); + ComponentMask component_mask(dim, default_component_mask); + + if (comp.find('x') != std::string::npos) + component_mask.set(0,true); + if (comp.find('y') != std::string::npos) + component_mask.set(1,true); + if (dim == 3 && comp.find('z') != std::string::npos) + component_mask.set(2,true); + + component_masks.push_back(component_mask); } } prm.leave_subsection(); // go through the list, create objects and let them parse // their own parameters - for (const auto &boundary_id : boundary_traction_indicators) + for (const auto &plugin_name: this->plugin_names) { - for (const auto &name : boundary_id.second.second) - { - boundary_traction_objects[boundary_id.first].push_back( - std::unique_ptr> (std::get(registered_plugins) - .create_plugin (name, - "Boundary traction::Model names"))); + // create boundary traction objects + this->plugin_objects.push_back(std::get(registered_plugins) + .create_plugin (plugin_name, + "Boundary traction::Model names")); - if (SimulatorAccess *sim = dynamic_cast*>(boundary_traction_objects[boundary_id.first].back().get())) - sim->initialize_simulator (this->get_simulator()); + if (SimulatorAccess *sim = dynamic_cast*>(this->plugin_objects.back().get())) + sim->initialize_simulator (this->get_simulator()); - boundary_traction_objects[boundary_id.first].back()->parse_parameters (prm); - boundary_traction_objects[boundary_id.first].back()->initialize (); - } + this->plugin_objects.back()->parse_parameters (prm); + this->plugin_objects.back()->initialize (); } } diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 83f5661a68a..eb949ac4612 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -939,9 +939,10 @@ namespace aspect const typename DoFHandler::face_iterator face = scratch.cell->face(scratch.face_number); - if (this->get_boundary_traction_manager().get_active_boundary_traction_names().find (face->boundary_id()) - != - this->get_boundary_traction_manager().get_active_boundary_traction_names().end()) + const auto &traction_bis = this->get_boundary_traction_manager().get_active_plugin_boundary_indicators(); + + if (std::find(traction_bis.begin(), traction_bis.end(), face->boundary_id()) + != traction_bis.end()) { for (unsigned int q=0; qstokes_system_on_boundary_face.push_back( std::make_unique>()); @@ -767,7 +767,7 @@ namespace aspect = ( // see if we need to assemble traction boundary conditions. // only if so do we actually need to have an FEFaceValues object - boundary_traction_manager.get_active_boundary_traction_names ().size() > 0 + !boundary_traction_manager.get_active_plugins().empty() ? update_values | update_quadrature_points | diff --git a/source/simulator/melt.cc b/source/simulator/melt.cc index c5110d07924..e478f3d7a5c 100644 --- a/source/simulator/melt.cc +++ b/source/simulator/melt.cc @@ -1034,11 +1034,10 @@ namespace aspect Assert(face_no != numbers::invalid_unsigned_int,ExcInternalError()); const typename DoFHandler::face_iterator face = cell->face(face_no); + const auto &traction_bis = this->get_boundary_traction_manager().get_active_plugin_boundary_indicators(); - if (this->get_boundary_traction_manager().get_active_boundary_traction_names() - .find (face->boundary_id()) - != - this->get_boundary_traction_manager().get_active_boundary_traction_names().end()) + if (std::find(traction_bis.begin(), traction_bis.end(), face->boundary_id()) + != traction_bis.end()) { scratch.face_finite_element_values.reinit (cell, face_no); @@ -1690,7 +1689,7 @@ namespace aspect std::make_unique>()); // add the terms for traction boundary conditions - if (!this->get_boundary_traction_manager().get_active_boundary_traction_names().empty()) + if (!this->get_boundary_traction_manager().get_active_plugins().empty()) { assemblers.stokes_system_on_boundary_face.push_back( std::make_unique> ()); diff --git a/source/simulator/newton.cc b/source/simulator/newton.cc index dbe07857332..e22adb8d381 100644 --- a/source/simulator/newton.cc +++ b/source/simulator/newton.cc @@ -97,7 +97,7 @@ namespace aspect " defined that handles this formulation.")); // add the terms for traction boundary conditions - if (!this->get_boundary_traction_manager().get_active_boundary_traction_names().empty()) + if (!this->get_boundary_traction_manager().get_active_plugins().empty()) { assemblers.stokes_system_on_boundary_face.push_back( std::make_unique>()); From ba097d3da4b03bd5040c2ebb414f6172a370a65f Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Tue, 22 Oct 2024 13:19:55 +0200 Subject: [PATCH 2/7] convert velocity boundary manager --- include/aspect/boundary_traction/interface.h | 12 + include/aspect/boundary_velocity/interface.h | 101 ++++++-- include/aspect/parameters.h | 8 - source/boundary_traction/interface.cc | 26 +- source/boundary_velocity/ascii_data.cc | 15 +- source/boundary_velocity/interface.cc | 116 ++++++--- source/mesh_deformation/free_surface.cc | 17 +- source/simulator/assembly.cc | 2 +- source/simulator/core.cc | 99 ++++--- source/simulator/helper_functions.cc | 255 ++++++++----------- source/simulator/solver_schemes.cc | 7 +- source/simulator/stokes_matrix_free.cc | 36 +-- tests/ascii_boundary_member.cc | 21 +- 13 files changed, 375 insertions(+), 340 deletions(-) diff --git a/include/aspect/boundary_traction/interface.h b/include/aspect/boundary_traction/interface.h index 85a8cdd13d6..a1235e0b073 100644 --- a/include/aspect/boundary_traction/interface.h +++ b/include/aspect/boundary_traction/interface.h @@ -119,6 +119,13 @@ namespace aspect const std::map>>> & get_active_boundary_traction_conditions () const; + /** + * Return a set of boundary indicators for which boundary + * tractions are prescribed. + */ + const std::set & + get_prescribed_boundary_traction_indicators () const; + /** * Return a list of boundary indicators that indicate for * each active plugin which boundary id @@ -220,6 +227,11 @@ namespace aspect */ std::vector component_masks; + /** + * A set of boundary indicators, on which tractions are prescribed. + */ + std::set prescribed_traction_boundary_indicators; + /** * A list of boundary traction objects that have been requested in the * parameter file. diff --git a/include/aspect/boundary_velocity/interface.h b/include/aspect/boundary_velocity/interface.h index 9527c9bf453..8f1ccf69fb5 100644 --- a/include/aspect/boundary_velocity/interface.h +++ b/include/aspect/boundary_velocity/interface.h @@ -79,28 +79,9 @@ namespace aspect * @ingroup BoundaryVelocities */ template - class Manager : public SimulatorAccess + class Manager : public Plugins::ManagerBase>, public SimulatorAccess { public: - /** - * Destructor. Made virtual since this class has virtual member - * functions. - */ - ~Manager () override; - - /** - * A function that is called at the beginning of each time step and - * calls the corresponding functions of all created plugins. - * - * The point of this function is to allow complex boundary velocity - * models to do an initialization step once at the beginning of each - * time step. An example would be a model that needs to call an - * external program to compute the velocity change at a boundary. - */ - virtual - void - update (); - /** * A function that calls the boundary_velocity functions of all the * individual boundary velocity objects and uses the stored operators @@ -146,7 +127,11 @@ namespace aspect * If there are no prescribed boundary velocity plugins * for a particular boundary, this boundary identifier will not appear * in the map. + * + * @deprecated: This function will be removed. Use the function + * get_active_plugin_names() of the base class ManagerBase instead. */ + DEAL_II_DEPRECATED const std::map>> & get_active_boundary_velocity_names () const; @@ -158,10 +143,39 @@ namespace aspect * boundary models for this boundary. If there are no prescribed * boundary velocity plugins for a particular boundary this boundary * identifier will not appear in the map. + * + * @deprecated: This function has been removed. Use the function + * get_active_plugins() of the base class ManagerBase instead. */ + DEAL_II_DEPRECATED const std::map>>> & get_active_boundary_velocity_conditions () const; + /** + * Return a list of boundary indicators that indicate for + * each active plugin which boundary id + * it is responsible for. The list of active plugins can be + * requested by calling get_active_plugins(). + */ + const std::vector & + get_active_plugin_boundary_indicators() const; + + /** + * Return a list of component masks that indicate for + * each active plugin which components it is responsible for. + * The list of plugin objects can be + * requested by calling get_active_plugins(). + */ + const std::vector & + get_active_plugin_component_masks() const; + + /** + * Return a set of boundary indicators for which boundary + * velocities are prescribed. + */ + const std::set & + get_prescribed_boundary_velocity_indicators () const; + /** * Return a list of boundary ids on which the velocity is prescribed * to be zero (no-slip). @@ -190,7 +204,7 @@ namespace aspect * then let these objects read their parameters as well. */ void - parse_parameters (ParameterHandler &prm); + parse_parameters (ParameterHandler &prm) override; /** * Go through the list of all boundary velocity models that have been selected @@ -200,9 +214,15 @@ namespace aspect * * This function can only be called if the given template type (the first template * argument) is a class derived from the Interface class in this namespace. + * + * @deprecated Instead of this function, use the + * Plugins::ManagerBase::has_matching_active_plugin() and + * Plugins::ManagerBase::get_matching_active_plugin() functions of the base + * class of the current class. */ template ,BoundaryVelocityType>::value>> + DEAL_II_DEPRECATED bool has_matching_boundary_velocity_model () const; @@ -216,9 +236,15 @@ namespace aspect * * This function can only be called if the given template type (the first template * argument) is a class derived from the Interface class in this namespace. + * + * @deprecated Instead of this function, use the + * Plugins::ManagerBase::has_matching_active_plugin() and + * Plugins::ManagerBase::get_matching_active_plugin() functions of the base + * class of the current class. */ template ,BoundaryVelocityType>::value>> + DEAL_II_DEPRECATED const BoundaryVelocityType & get_matching_boundary_velocity_model () const; @@ -245,9 +271,32 @@ namespace aspect << arg1 << "> among the names of registered boundary velocity objects."); private: + /** + * A list of boundary indicators that indicate for + * each plugin in the list of plugin_objects which boundary id + * it is responsible for. By default each plugin + * is active for all boundaries, but this list + * can be modified by derived classes to limit the application + * of plugins to specific boundaries. + */ + std::vector boundary_indicators; + + /** + * A list of boundary indicators that indicate for + * each plugin in the list of plugin_objects which components + * it is responsible for. By default each plugin + * is active for all components, but this list + * can be modified by derived classes to limit the application + * of plugins to specific boundaries. + */ + std::vector component_masks; + /** * A list of boundary velocity objects that have been requested in the * parameter file. + * + * @deprecated: This variable is no longer used, but needed to issue a proper + * error message in the function get_active_boundary_velocity_conditions(). */ std::map>>> boundary_velocity_objects; @@ -258,9 +307,19 @@ namespace aspect * mapped to one of the plugins of velocity boundary conditions (e.g. * "function"). If the components string is empty, it is assumed the * plugins are used for all components. + * + * @deprecated: Remove this variable when the deprecated functions + * get_active_boundary_velocity_names and + * get_active_boundary_velocity_conditions are removed. Use the base class + * variable plugin_names instead. */ std::map>> boundary_velocity_indicators; + /** + * A set of boundary indicators, on which velocities are prescribed. + */ + std::set prescribed_velocity_boundary_indicators; + /** * A set of boundary indicators, on which velocities are prescribed to * zero (no-slip). diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h index 7fdfb7aa265..4c29d9e692c 100644 --- a/include/aspect/parameters.h +++ b/include/aspect/parameters.h @@ -637,14 +637,6 @@ namespace aspect bool enable_additional_stokes_rhs; bool enable_prescribed_dilation; - /** - * Map from boundary id to a pair "components", "traction boundary type", - * where components is of the format "[x][y][z]" and the traction type is - * mapped to one of the plugins of traction boundary conditions (e.g. - * "function") - */ - std::map> prescribed_traction_boundary_indicators; - /** * A set of boundary ids on which the boundary_heat_flux objects * will be applied. diff --git a/source/boundary_traction/interface.cc b/source/boundary_traction/interface.cc index 1917009295d..c234437feb5 100644 --- a/source/boundary_traction/interface.cc +++ b/source/boundary_traction/interface.cc @@ -113,6 +113,15 @@ namespace aspect + template + const std::set & + Manager::get_prescribed_boundary_traction_indicators () const + { + return prescribed_traction_boundary_indicators; + } + + + template const std::vector & Manager::get_active_plugin_boundary_indicators() const @@ -273,16 +282,17 @@ namespace aspect this->plugin_names.push_back(value); boundary_indicators.push_back(boundary_id); + prescribed_traction_boundary_indicators.insert(boundary_id); - const bool default_component_mask = comp.empty(); - ComponentMask component_mask(dim, default_component_mask); + ComponentMask component_mask(this->introspection().n_components, + false); - if (comp.find('x') != std::string::npos) - component_mask.set(0,true); - if (comp.find('y') != std::string::npos) - component_mask.set(1,true); - if (dim == 3 && comp.find('z') != std::string::npos) - component_mask.set(2,true); + if (comp.empty() || comp.find('x') != std::string::npos) + component_mask.set(this->introspection().component_indices.velocities[0],true); + if (comp.empty() || comp.find('y') != std::string::npos) + component_mask.set(this->introspection().component_indices.velocities[1],true); + if (dim == 3 && (comp.empty() || comp.find('z') != std::string::npos)) + component_mask.set(this->introspection().component_indices.velocities[2],true); component_masks.push_back(component_mask); } diff --git a/source/boundary_velocity/ascii_data.cc b/source/boundary_velocity/ascii_data.cc index 5e0843075e1..faab1e46066 100644 --- a/source/boundary_velocity/ascii_data.cc +++ b/source/boundary_velocity/ascii_data.cc @@ -38,14 +38,17 @@ namespace aspect void AsciiData::initialize () { - for (const auto &p : this->get_boundary_velocity_manager().get_active_boundary_velocity_conditions()) + unsigned int i=0; + for (const auto &plugin : this->get_boundary_velocity_manager().get_active_plugins()) { - for (const auto &plugin : p.second) - if (plugin.get() == this) - boundary_ids.insert(p.first); + if (plugin.get() == this) + boundary_ids.insert(this->get_boundary_velocity_manager().get_active_plugin_boundary_indicators()[i]); + + ++i; } - AssertThrow(*(boundary_ids.begin()) != numbers::invalid_boundary_id, - ExcMessage("Did not find the boundary indicator for the prescribed data plugin.")); + + AssertThrow(boundary_ids.empty() == false, + ExcMessage("Did not find the boundary indicator for the velocity ascii data plugin.")); Utilities::AsciiDataBoundary::initialize(boundary_ids, dim); diff --git a/source/boundary_velocity/interface.cc b/source/boundary_velocity/interface.cc index a80a9245ad4..12b9ad48ba8 100644 --- a/source/boundary_velocity/interface.cc +++ b/source/boundary_velocity/interface.cc @@ -36,24 +36,6 @@ namespace aspect // ------------------------------ Manager ----------------------------- // -------------------------------- Deal with registering boundary_velocity models and automating // -------------------------------- their setup and selection at run time - - template - Manager::~Manager() - = default; - - - - template - void - Manager::update () - { - for (const auto &boundary : boundary_velocity_objects) - for (const auto &p : boundary.second) - p->update(); - } - - - namespace { std::tuple @@ -84,26 +66,38 @@ namespace aspect Manager::boundary_velocity (const types::boundary_id boundary_indicator, const Point &position) const { - typename std::map>>>::const_iterator boundary_plugins = - boundary_velocity_objects.find(boundary_indicator); + Tensor<1,dim> velocity; - Assert(boundary_plugins != boundary_velocity_objects.end(), + bool found_plugin = false; + unsigned int i=0; + for (const auto &plugin: this->plugin_objects) + { + if (boundary_indicators[i] == boundary_indicator) + { + found_plugin = true; + const Tensor<1,dim> plugin_velocity = plugin->boundary_velocity(boundary_indicator, + position); + for (unsigned int d=0; d velocity = Tensor<1,dim>(); - - for (const auto &plugin : boundary_plugins->second) - velocity += plugin->boundary_velocity(boundary_indicator, - position); - return velocity; } template + DEAL_II_DEPRECATED const std::map>> & Manager::get_active_boundary_velocity_names () const { @@ -113,14 +107,45 @@ namespace aspect template + DEAL_II_DEPRECATED const std::map>>> & Manager::get_active_boundary_velocity_conditions () const { + AssertThrow(false, ExcMessage("This function has been removed. Use the function " + "get_active_plugins() of the base class ManagerBase " + "instead.")); return boundary_velocity_objects; } + template + const std::vector & + Manager::get_active_plugin_boundary_indicators() const + { + return boundary_indicators; + } + + + + template + const std::vector & + Manager::get_active_plugin_component_masks() const + { + return component_masks; + } + + + + template + const std::set & + Manager::get_prescribed_boundary_velocity_indicators () const + { + return prescribed_velocity_boundary_indicators; + } + + + template const std::set & Manager::get_zero_boundary_velocity_indicators () const @@ -349,27 +374,40 @@ namespace aspect { boundary_velocity_indicators[boundary_id] = std::make_pair(comp,std::vector(1,value)); } + + this->plugin_names.push_back(value); + boundary_indicators.push_back(boundary_id); + prescribed_velocity_boundary_indicators.insert(boundary_id); + + ComponentMask component_mask(this->introspection().n_components, + false); + + if (comp.empty() || comp.find('x') != std::string::npos) + component_mask.set(this->introspection().component_indices.velocities[0],true); + if (comp.empty() || comp.find('y') != std::string::npos) + component_mask.set(this->introspection().component_indices.velocities[1],true); + if (dim == 3 && (comp.empty() || comp.find('z') != std::string::npos)) + component_mask.set(this->introspection().component_indices.velocities[2],true); + + component_masks.push_back(component_mask); } } prm.leave_subsection(); // go through the list, create objects and let them parse // their own parameters - for (const auto &boundary_id : boundary_velocity_indicators) + for (const auto &plugin_name: this->plugin_names) { - for (const auto &name : boundary_id.second.second) - { - boundary_velocity_objects[boundary_id.first].push_back( - std::unique_ptr> (std::get(registered_plugins) - .create_plugin (name, - "Boundary velocity::Model names"))); + // create boundary traction objects + this->plugin_objects.push_back(std::get(registered_plugins) + .create_plugin (plugin_name, + "Boundary velocity::Model names")); - if (SimulatorAccess *sim = dynamic_cast*>(boundary_velocity_objects[boundary_id.first].back().get())) - sim->initialize_simulator (this->get_simulator()); + if (SimulatorAccess *sim = dynamic_cast*>(this->plugin_objects.back().get())) + sim->initialize_simulator (this->get_simulator()); - boundary_velocity_objects[boundary_id.first].back()->parse_parameters (prm); - boundary_velocity_objects[boundary_id.first].back()->initialize (); - } + this->plugin_objects.back()->parse_parameters (prm); + this->plugin_objects.back()->initialize (); } } diff --git a/source/mesh_deformation/free_surface.cc b/source/mesh_deformation/free_surface.cc index 8d67a0f0557..b1ae199ac5e 100644 --- a/source/mesh_deformation/free_surface.cc +++ b/source/mesh_deformation/free_surface.cc @@ -50,20 +50,17 @@ namespace aspect std::set velocity_boundary_indicators = this->get_boundary_velocity_manager().get_zero_boundary_velocity_indicators(); // Get the tangential velocity boundary indicators - const std::set tmp_tangential_vel_boundary_indicators = this->get_boundary_velocity_manager().get_tangential_boundary_velocity_indicators(); - velocity_boundary_indicators.insert(tmp_tangential_vel_boundary_indicators.begin(), - tmp_tangential_vel_boundary_indicators.end()); + const auto &tangential_boundary_indicators = this->get_boundary_velocity_manager().get_tangential_boundary_velocity_indicators(); + velocity_boundary_indicators.insert(tangential_boundary_indicators.begin(), + tangential_boundary_indicators.end()); // Get the active velocity boundary indicators - const std::map>> - tmp_active_vel_boundary_indicators = this->get_boundary_velocity_manager().get_active_boundary_velocity_names(); - - for (const auto &p : tmp_active_vel_boundary_indicators) - velocity_boundary_indicators.insert(p.first); + const auto &prescribed_boundary_indicators = this->get_boundary_velocity_manager().get_prescribed_boundary_velocity_indicators(); + velocity_boundary_indicators.insert(prescribed_boundary_indicators.begin(), + prescribed_boundary_indicators.end()); // Get the mesh deformation boundary indicators - const std::set tmp_mesh_deformation_boundary_indicators = this->get_mesh_deformation_boundary_indicators(); - for (const auto &p : tmp_mesh_deformation_boundary_indicators) + for (const auto &p : this->get_mesh_deformation_boundary_indicators()) AssertThrow(velocity_boundary_indicators.find(p) == velocity_boundary_indicators.end(), ExcMessage("The free surface mesh deformation plugin cannot be used with the current velocity boundary conditions")); } diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc index bb54b960cb6..289c6f1b5c2 100644 --- a/source/simulator/assembly.cc +++ b/source/simulator/assembly.cc @@ -743,7 +743,7 @@ namespace aspect // we will update the right-hand side with boundary information in // StokesMatrixFreeHandler::correct_stokes_rhs(). if (!stokes_matrix_free) - Assert(rebuild_stokes_matrix || boundary_velocity_manager.get_active_boundary_velocity_conditions().size()==0, + Assert(rebuild_stokes_matrix || boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty(), ExcInternalError("If we have inhomogeneous constraints, we must re-assemble the system matrix.")); system_rhs = 0; diff --git a/source/simulator/core.cc b/source/simulator/core.cc index e8173dae562..515478a5f88 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -505,8 +505,8 @@ namespace aspect std::set open_velocity_boundary_indicators = geometry_model->get_used_boundary_indicators(); - for (const auto &p : boundary_velocity_manager.get_active_boundary_velocity_names()) - open_velocity_boundary_indicators.erase (p.first); + for (const auto p : boundary_velocity_manager.get_prescribed_boundary_velocity_indicators()) + open_velocity_boundary_indicators.erase (p); for (const auto p : boundary_velocity_manager.get_zero_boundary_velocity_indicators()) open_velocity_boundary_indicators.erase (p); for (const auto p : boundary_velocity_manager.get_tangential_boundary_velocity_indicators()) @@ -1375,67 +1375,54 @@ namespace aspect // set the current time and do the interpolation // for the prescribed velocity fields boundary_velocity_manager.update(); - for (const auto &p : boundary_velocity_manager.get_active_boundary_velocity_names()) + + // translate from component mask per plugin to component mask per boundary id + std::map boundary_component_masks; + for (unsigned int i=0; isecond == component_mask, + ExcMessage("Boundary indicator <" + + + Utilities::int_to_string(boundary_id) + + + "> with symbolic name <" + + + geometry_model->translate_id_to_symbol_name (boundary_id) + + + "> is listed as having different component masks " + "for different boundary velocity plugins. This is not allowed.")); + } + else + boundary_component_masks[boundary_id] = component_mask; + } + + // put boundary conditions into constraints object for each boundary + for (const auto boundary_id: boundary_velocity_manager.get_prescribed_boundary_velocity_indicators()) { Utilities::VectorFunctionFromVelocityFunctionObject vel (introspection.n_components, [&] (const dealii::Point &x) -> Tensor<1,dim> { - return boundary_velocity_manager.boundary_velocity(p.first, x); - }); - - // here we create a mask for interpolate_boundary_values out of the 'selector' - std::vector mask(introspection.component_masks.velocities.size(), false); - const std::string &comp = p.second.first; + if (!assemble_newton_stokes_system || (assemble_newton_stokes_system && nonlinear_iteration == 0)) + return boundary_velocity_manager.boundary_velocity(boundary_id, x); + else + return Tensor<1,dim>(); - if (comp.length()>0) - { - for (const char direction : comp) - { - switch (direction) - { - case 'x': - mask[introspection.component_indices.velocities[0]] = true; - break; - case 'y': - mask[introspection.component_indices.velocities[1]] = true; - break; - case 'z': - // we must be in 3d, or 'z' should never have gotten through - Assert (dim==3, ExcInternalError()); - if (dim==3) - mask[introspection.component_indices.velocities[dim-1]] = true; - break; - default: - Assert (false, ExcInternalError()); - } - } - } - else - { - // no mask given -- take all velocities - for (unsigned int i=0; i(); + }); - if (!assemble_newton_stokes_system || (assemble_newton_stokes_system && nonlinear_iteration == 0)) - { - VectorTools::interpolate_boundary_values (*mapping, - dof_handler, - p.first, - vel, - constraints, - ComponentMask(mask)); - } - else - { - VectorTools::interpolate_boundary_values (*mapping, - dof_handler, - p.first, - Functions::ZeroFunction(introspection.n_components), - constraints, - ComponentMask(mask)); - } + VectorTools::interpolate_boundary_values (*mapping, + dof_handler, + boundary_id, + vel, + constraints, + boundary_component_masks.at(boundary_id)); } } diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc index 8a5a3742909..d54b5ee18a2 100644 --- a/source/simulator/helper_functions.cc +++ b/source/simulator/helper_functions.cc @@ -2290,7 +2290,7 @@ namespace aspect boundary_velocity_manager.get_zero_boundary_velocity_indicators(); const auto &prescribed_velocity_boundaries = - boundary_velocity_manager.get_active_boundary_velocity_conditions(); + boundary_velocity_manager.get_prescribed_boundary_velocity_indicators(); // Loop over all of the boundary faces, ... for (const auto &cell : dof_handler.active_cell_iterators()) @@ -2380,104 +2380,50 @@ namespace aspect void Simulator::check_consistency_of_boundary_conditions() const { - // make sure velocity and traction boundary indicators don't appear in multiple lists - std::set boundary_indicator_lists[6] - = { boundary_velocity_manager.get_zero_boundary_velocity_indicators(), - boundary_velocity_manager.get_tangential_boundary_velocity_indicators(), - std::set() // to be prescribed velocity and traction boundary indicators - }; - - // sets of the boundary indicators only (no selectors and values) - std::set velocity_bi; - std::set traction_bi; - - for (const auto &p : boundary_velocity_manager.get_active_boundary_velocity_names()) - velocity_bi.insert(p.first); - - for (const auto &r : parameters.prescribed_traction_boundary_indicators) - traction_bi.insert(r.first); - - // are there any indicators that occur in both the prescribed velocity and traction list? - std::set intersection; - std::set_intersection (velocity_bi.begin(), - velocity_bi.end(), - traction_bi.begin(), - traction_bi.end(), - std::inserter(intersection, intersection.end())); - - // if so, do they have different selectors? - if (!intersection.empty()) - { - for (const auto it : intersection) - { - const std::map>>::const_iterator - boundary_velocity_names = boundary_velocity_manager.get_active_boundary_velocity_names().find(it); - Assert(boundary_velocity_names != boundary_velocity_manager.get_active_boundary_velocity_names().end(), - ExcInternalError()); - - std::set velocity_selector; - std::set traction_selector; - - for (const auto it_selector : boundary_velocity_names->second.first) - velocity_selector.insert(it_selector); - - for (std::string::const_iterator - it_selector = parameters.prescribed_traction_boundary_indicators.find(it)->second.first.begin(); - it_selector != parameters.prescribed_traction_boundary_indicators.find(it)->second.first.end(); - ++it_selector) - traction_selector.insert(*it_selector); - - // if there are no selectors specified, throw exception - AssertThrow(!velocity_selector.empty() || !traction_selector.empty(), - ExcMessage ("Boundary indicator <" - + - Utilities::int_to_string(it) - + - "> with symbolic name <" - + - geometry_model->translate_id_to_symbol_name (it) - + - "> is listed as having both " - "velocity and traction boundary conditions in the input file.")); - - std::set intersection_selector; - std::set_intersection (velocity_selector.begin(), - velocity_selector.end(), - traction_selector.begin(), - traction_selector.end(), - std::inserter(intersection_selector, intersection_selector.end())); - - // if the same selectors are specified, throw exception - AssertThrow(intersection_selector.empty(), - ExcMessage ("Selectors of boundary indicator <" - + - Utilities::int_to_string(it) - + - "> with symbolic name <" - + - geometry_model->translate_id_to_symbol_name (it) - + - "> are listed as having both " - "velocity and traction boundary conditions in the input file.")); - } - } - - // remove correct boundary indicators that occur in both the velocity and the traction set - // but have different selectors - std::set union_set; - std::set_union (velocity_bi.begin(), - velocity_bi.end(), - traction_bi.begin(), - traction_bi.end(), - std::inserter(union_set, union_set.end())); - - // assign the prescribed boundary indicator list to the boundary_indicator_lists - boundary_indicator_lists[3] = union_set; + // a container for the indicators of all boundary conditions + std::vector> boundary_indicator_lists; + boundary_indicator_lists.emplace_back(boundary_velocity_manager.get_zero_boundary_velocity_indicators()); + boundary_indicator_lists.emplace_back(boundary_velocity_manager.get_tangential_boundary_velocity_indicators()); + boundary_indicator_lists.emplace_back(boundary_velocity_manager.get_prescribed_boundary_velocity_indicators()); + boundary_indicator_lists.emplace_back(boundary_traction_manager.get_prescribed_boundary_traction_indicators()); + + // Make sure that each combination of boundary velocity and boundary traction plugin + // either refers to different boundary indicators or to different components + const auto &velocity_boundary_ids = boundary_velocity_manager.get_active_plugin_boundary_indicators(); + const auto &traction_boundary_ids = boundary_traction_manager.get_active_plugin_boundary_indicators(); + + for (unsigned int i=0; i with symbolic name <" + + + geometry_model->translate_id_to_symbol_name (velocity_boundary_ids[i]) + + + "> is listed as having both " + "velocity and traction boundary conditions in the input file.")); + + // we have ensured the prescribed velocity and prescribed traction boundary conditions + // are compatible. In order to check them against the other boundary conditions, we + // need to remove the boundary indicator from one of the lists to make sure it only + // appears in one of them. We choose to remove the boundary indicator from the traction list. + boundary_indicator_lists[3].erase(traction_boundary_ids[j]); + } + } - // for each combination of boundary indicator lists, make sure that the + // for each combination of velocity boundary indicator lists, make sure that the // intersection is empty - for (unsigned int i=0; i intersection; std::set_intersection (boundary_indicator_lists[i].begin(), @@ -2502,20 +2448,28 @@ namespace aspect // make sure temperature and heat flux boundary indicators don't appear in multiple lists // this is easier than for the velocity/traction, as there are no selectors - std::set temperature_bi = boundary_temperature_manager.get_fixed_temperature_boundary_indicators(); - std::set heat_flux_bi = parameters.fixed_heat_flux_boundary_indicators; + boundary_indicator_lists.emplace_back(boundary_temperature_manager.get_fixed_temperature_boundary_indicators()); + boundary_indicator_lists.emplace_back(parameters.fixed_heat_flux_boundary_indicators); // are there any indicators that occur in both the prescribed temperature and heat flux list? std::set T_intersection; - std::set_intersection (temperature_bi.begin(), - temperature_bi.end(), - heat_flux_bi.begin(), - heat_flux_bi.end(), + std::set_intersection (boundary_temperature_manager.get_fixed_temperature_boundary_indicators().begin(), + boundary_temperature_manager.get_fixed_temperature_boundary_indicators().end(), + parameters.fixed_heat_flux_boundary_indicators.begin(), + parameters.fixed_heat_flux_boundary_indicators.end(), std::inserter(T_intersection, T_intersection.end())); - AssertThrow(T_intersection.empty(), - ExcMessage ("There is a boundary indicator listed as having both " - "temperature and heat flux boundary conditions in the input file.")); + AssertThrow (T_intersection.empty(), + ExcMessage ("Boundary indicator <" + + + Utilities::int_to_string(*T_intersection.begin()) + + + "> with symbolic name <" + + + geometry_model->translate_id_to_symbol_name (*T_intersection.begin()) + + + "> is listed as having more " + "than one type of temperature or heat flux boundary condition in the input file.")); // Check that the periodic boundaries do not have other boundary conditions set using periodic_boundary_set @@ -2524,56 +2478,55 @@ namespace aspect periodic_boundary_set pbs = geometry_model->get_periodic_boundary_pairs(); for (const auto &pb : pbs) - { - // Throw error if we are trying to use the same boundary for more than one boundary condition - AssertThrow( is_element( pb.first.first, boundary_temperature_manager.get_fixed_temperature_boundary_indicators() ) == false && - is_element( pb.first.second, boundary_temperature_manager.get_fixed_temperature_boundary_indicators() ) == false && - is_element( pb.first.first, boundary_composition_manager.get_fixed_composition_boundary_indicators() ) == false && - is_element( pb.first.second, boundary_composition_manager.get_fixed_composition_boundary_indicators() ) == false && - is_element( pb.first.first, boundary_indicator_lists[0] ) == false && // zero velocity - is_element( pb.first.second, boundary_indicator_lists[0] ) == false && // zero velocity - is_element( pb.first.first, boundary_indicator_lists[1] ) == false && // tangential velocity - is_element( pb.first.second, boundary_indicator_lists[1] ) == false && // tangential velocity - is_element( pb.first.first, boundary_indicator_lists[3] ) == false && // prescribed traction or velocity - is_element( pb.first.second, boundary_indicator_lists[3] ) == false, // prescribed traction or velocity - ExcMessage("Periodic boundaries must not have boundary conditions set.")); - } + for (const auto &boundary_indicators: boundary_indicator_lists) + { + AssertThrow(is_element(pb.first.first, boundary_indicators) == false, + ExcMessage ("Boundary indicator <" + + + Utilities::int_to_string(pb.first.first) + + + "> with symbolic name <" + + + geometry_model->translate_id_to_symbol_name (pb.first.first) + + + "> is listed as having a periodic boundary condition " + "in the input file, but also has another type of boundary condition. " + "Periodic boundaries cannot have other boundary conditions.")); + + AssertThrow(is_element(pb.first.second, boundary_indicators) == false, + ExcMessage ("Boundary indicator <" + + + Utilities::int_to_string(pb.first.second) + + + "> with symbolic name <" + + + geometry_model->translate_id_to_symbol_name (pb.first.second) + + + "> is listed as having a periodic boundary condition " + "in the input file, but also has another type of boundary condition. " + "Periodic boundaries cannot have other boundary conditions.")); + } const std::set all_boundary_indicators = geometry_model->get_used_boundary_indicators(); - if (parameters.nonlinear_solver != NonlinearSolver::single_Advection_no_Stokes) - { - // next make sure that all listed indicators are actually used by - // this geometry - for (const auto &list : boundary_indicator_lists) - for (const auto &p : list) - AssertThrow (all_boundary_indicators.find (p) - != all_boundary_indicators.end(), - ExcMessage ("One of the boundary indicators listed in the input file " - "is not used by the geometry model.")); - } - else + + // next make sure that all listed indicators are actually used by + // this geometry + for (const auto &list : boundary_indicator_lists) + for (const auto &p : list) + AssertThrow (all_boundary_indicators.find (p) + != all_boundary_indicators.end(), + ExcMessage ("One of the boundary indicators listed in the input file " + "is not used by the geometry model.")); + + if (parameters.nonlinear_solver == NonlinearSolver::single_Advection_no_Stokes) { - // next make sure that there are no listed indicators - for (const auto &list : boundary_indicator_lists) - AssertThrow (list.empty(), + // make sure that there are no listed velocity boundary conditions + for (unsigned int i=0; i<4; ++i) + AssertThrow (boundary_indicator_lists[i].empty(), ExcMessage ("With the solver scheme `single Advection, no Stokes', " "one cannot set boundary conditions for velocity.")); } - - - // now do the same for the fixed temperature indicators and the - // compositional indicators - for (const auto p : boundary_temperature_manager.get_fixed_temperature_boundary_indicators()) - AssertThrow (all_boundary_indicators.find (p) - != all_boundary_indicators.end(), - ExcMessage ("One of the fixed boundary temperature indicators listed in the input file " - "is not used by the geometry model.")); - for (const auto p : boundary_composition_manager.get_fixed_composition_boundary_indicators()) - AssertThrow (all_boundary_indicators.find (p) - != all_boundary_indicators.end(), - ExcMessage ("One of the fixed boundary composition indicators listed in the input file " - "is not used by the geometry model.")); } @@ -2598,8 +2551,8 @@ namespace aspect // rebuild the whole system to compute the rhs. assemble_newton_stokes_system = true; rebuild_stokes_preconditioner = false; - rebuild_stokes_matrix = boundary_velocity_manager.get_active_boundary_velocity_conditions().size()!=0; - assemble_newton_stokes_matrix = boundary_velocity_manager.get_active_boundary_velocity_conditions().size()!=0; + rebuild_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty(); + assemble_newton_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty(); compute_current_constraints (); diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc index 4c4f71ba4a9..e21aaa61016 100644 --- a/source/simulator/solver_schemes.cc +++ b/source/simulator/solver_schemes.cc @@ -370,7 +370,7 @@ namespace aspect // don't need to force assembly of the matrix. if (stokes_matrix_depends_on_solution() || - (boundary_velocity_manager.get_active_boundary_velocity_conditions().size() > 0) + (boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().size() > 0) || parameters.mesh_deformation_enabled ) rebuild_stokes_matrix = rebuild_stokes_preconditioner = true; @@ -579,7 +579,7 @@ namespace aspect // don't need to force assembly of the matrix. if (stokes_matrix_depends_on_solution() || - (nonlinear_iteration == 0 && boundary_velocity_manager.get_active_boundary_velocity_conditions().size() > 0)) + (nonlinear_iteration == 0 && boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().size() > 0)) rebuild_stokes_matrix = rebuild_stokes_preconditioner = assemble_newton_stokes_matrix = true; else if (parameters.enable_prescribed_dilation) // The dilation requires the Stokes matrix (which is on the rhs @@ -680,8 +680,7 @@ namespace aspect // Rebuild the rhs to determine the new residual. assemble_newton_stokes_matrix = rebuild_stokes_preconditioner = false; - rebuild_stokes_matrix = (boundary_velocity_manager.get_active_boundary_velocity_conditions().empty() - == false); + rebuild_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty(); assemble_stokes_system(); diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc index f867d462cb4..5a716de77ce 100644 --- a/source/simulator/stokes_matrix_free.cc +++ b/source/simulator/stokes_matrix_free.cc @@ -2487,35 +2487,19 @@ namespace aspect mg_constrained_dofs_A_block.initialize(dof_handler_v); std::set dirichlet_boundary = sim.boundary_velocity_manager.get_zero_boundary_velocity_indicators(); - for (const auto &it: sim.boundary_velocity_manager.get_active_boundary_velocity_names()) + for (unsigned int i=0; i0) + if (component_mask != ComponentMask(sim.introspection.n_components, false)) { - std::vector mask(fe_v.n_components(), false); - for (const auto &direction : component) - { - switch (direction) - { - case 'x': - mask[0] = true; - break; - case 'y': - mask[1] = true; - break; - case 'z': - // we must be in 3d, or 'z' should never have gotten through - Assert (dim==3, ExcInternalError()); - if (dim==3) - mask[2] = true; - break; - default: - Assert (false, ExcInternalError()); - } - } - mg_constrained_dofs_A_block.make_zero_boundary_constraints(dof_handler_v, {bdryid}, ComponentMask(mask)); + ComponentMask mask(fe_v.n_components(), false); + + for (unsigned int i=0; i::initialize () { - const std::map>>> & - bvs = this->get_boundary_velocity_manager().get_active_boundary_velocity_conditions(); - for (const auto &boundary : bvs) - for (const auto &p : boundary.second) - { - if (p.get() == this) - boundary_ids.insert(boundary.first); - } - AssertThrow(*(boundary_ids.begin()) != numbers::invalid_boundary_id, - ExcMessage("Did not find the boundary indicator for the prescribed data plugin.")); + unsigned int i=0; + for (const auto &plugin : this->get_boundary_velocity_manager().get_active_plugins()) + { + if (plugin.get() == this) + boundary_ids.insert(this->get_boundary_velocity_manager().get_active_plugin_boundary_indicators()[i]); + + ++i; + } + + AssertThrow(boundary_ids.empty() == false, + ExcMessage("Did not find the boundary indicator for the velocity ascii data plugin.")); member->initialize(boundary_ids, dim); From cd5d2ae95227a42adee55e2eb60553f6a02cfd1a Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 23 Oct 2024 10:50:37 +0200 Subject: [PATCH 3/7] Update tests --- tests/ascii_boundary_member/screen-output | 6 +++--- .../screen-output | 6 +++--- .../screen-output | 12 ++++++------ .../screen-output | 8 ++++---- .../screen-output | 8 ++++---- .../screen-output | 8 ++++---- .../screen-output | 12 ++++++------ .../screen-output | 10 +++++----- .../screen-output | 10 +++++----- .../screen-output | 10 +++++----- 10 files changed, 45 insertions(+), 45 deletions(-) diff --git a/tests/ascii_boundary_member/screen-output b/tests/ascii_boundary_member/screen-output index ee25bc51e2a..65e032a82ed 100644 --- a/tests/ascii_boundary_member/screen-output +++ b/tests/ascii_boundary_member/screen-output @@ -2,13 +2,13 @@ Loading shared library <./libascii_boundary_member.debug.so> - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_left.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_bottom.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_right.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_left.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_bottom.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_right.0.txt. Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_top.0.txt. diff --git a/tests/ascii_data_boundary_velocity_2d_box/screen-output b/tests/ascii_data_boundary_velocity_2d_box/screen-output index cabba6ace67..ac1fa9fc1fc 100644 --- a/tests/ascii_data_boundary_velocity_2d_box/screen-output +++ b/tests/ascii_data_boundary_velocity_2d_box/screen-output @@ -1,12 +1,12 @@ - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_left.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_bottom.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_right.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_left.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_bottom.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_right.0.txt. Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_top.0.txt. diff --git a/tests/ascii_data_boundary_velocity_2d_box_time/screen-output b/tests/ascii_data_boundary_velocity_2d_box_time/screen-output index f29e4ee7c68..3895a70a931 100644 --- a/tests/ascii_data_boundary_velocity_2d_box_time/screen-output +++ b/tests/ascii_data_boundary_velocity_2d_box_time/screen-output @@ -1,5 +1,11 @@ + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_top.0.txt. + + + Also loading next Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_top.1.txt. + + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_left.0.txt. @@ -21,12 +27,6 @@ that time-dependence ends at this timestep (i.e. the boundary condition will continue unchanged from the last known state into the future). - - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_top.0.txt. - - - Also loading next Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_2d_top.1.txt. - Number of active cells: 80 (on 3 levels) Number of degrees of freedom: 1,212 (738+105+369) diff --git a/tests/ascii_data_boundary_velocity_2d_chunk/screen-output b/tests/ascii_data_boundary_velocity_2d_chunk/screen-output index 4317af29f5b..d526dc586f6 100644 --- a/tests/ascii_data_boundary_velocity_2d_chunk/screen-output +++ b/tests/ascii_data_boundary_velocity_2d_chunk/screen-output @@ -1,15 +1,15 @@ - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/chunk_2d_bottom.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/chunk_2d_west.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/chunk_2d_top.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/chunk_2d_east.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/chunk_2d_west.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/chunk_2d_bottom.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/chunk_2d_east.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/chunk_2d_top.0.txt. Number of active cells: 64 (on 4 levels) Number of degrees of freedom: 948 (578+81+289) diff --git a/tests/ascii_data_boundary_velocity_2d_shell/screen-output b/tests/ascii_data_boundary_velocity_2d_shell/screen-output index 32529ac4498..121fad88cea 100644 --- a/tests/ascii_data_boundary_velocity_2d_shell/screen-output +++ b/tests/ascii_data_boundary_velocity_2d_shell/screen-output @@ -1,15 +1,15 @@ - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_bottom.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_left.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_top.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_right.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_left.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_bottom.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_right.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_top.0.txt. Number of active cells: 192 (on 4 levels) Number of degrees of freedom: 2,724 (1,666+225+833) diff --git a/tests/ascii_data_boundary_velocity_2d_shell_spherical/screen-output b/tests/ascii_data_boundary_velocity_2d_shell_spherical/screen-output index 8081720c373..e8778867c3a 100644 --- a/tests/ascii_data_boundary_velocity_2d_shell_spherical/screen-output +++ b/tests/ascii_data_boundary_velocity_2d_shell_spherical/screen-output @@ -1,15 +1,15 @@ - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_bottom_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_left_spherical.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_top_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_right_spherical.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_left_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_bottom_spherical.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_right_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_2d_top_spherical.0.txt. Number of active cells: 192 (on 4 levels) Number of degrees of freedom: 2,724 (1,666+225+833) diff --git a/tests/ascii_data_boundary_velocity_3d_box/screen-output b/tests/ascii_data_boundary_velocity_3d_box/screen-output index bca8872f2f0..ae77e226ff7 100644 --- a/tests/ascii_data_boundary_velocity_3d_box/screen-output +++ b/tests/ascii_data_boundary_velocity_3d_box/screen-output @@ -1,21 +1,21 @@ - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_left.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_bottom.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_right.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_left.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_front.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_right.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_back.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_top.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_bottom.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_front.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_top.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/box_3d_back.0.txt. Number of active cells: 200 (on 2 levels) Number of degrees of freedom: 9,183 (6,615+363+2,205) diff --git a/tests/ascii_data_boundary_velocity_3d_shell/screen-output b/tests/ascii_data_boundary_velocity_3d_shell/screen-output index 8cf41255551..bb70a37f1e0 100644 --- a/tests/ascii_data_boundary_velocity_3d_shell/screen-output +++ b/tests/ascii_data_boundary_velocity_3d_shell/screen-output @@ -1,18 +1,18 @@ - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_bottom.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_west.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_top.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_east.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_east.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_south.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_west.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_bottom.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_south.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_top.0.txt. Number of active cells: 24 (on 2 levels) Number of degrees of freedom: 1,277 (915+57+305) diff --git a/tests/ascii_data_boundary_velocity_3d_shell_spherical/screen-output b/tests/ascii_data_boundary_velocity_3d_shell_spherical/screen-output index 1ff28fafe99..ad4ea20d99b 100644 --- a/tests/ascii_data_boundary_velocity_3d_shell_spherical/screen-output +++ b/tests/ascii_data_boundary_velocity_3d_shell_spherical/screen-output @@ -1,18 +1,18 @@ - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_bottom_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_west_spherical.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_top_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_east_spherical.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_east_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_south_spherical.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_west_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_bottom_spherical.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_south_spherical.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_top_spherical.0.txt. Number of active cells: 24 (on 2 levels) Number of degrees of freedom: 1,277 (915+57+305) diff --git a/tests/ascii_data_boundary_velocity_residual_shell/screen-output b/tests/ascii_data_boundary_velocity_residual_shell/screen-output index 3334c2e411e..6521a88cdab 100644 --- a/tests/ascii_data_boundary_velocity_residual_shell/screen-output +++ b/tests/ascii_data_boundary_velocity_residual_shell/screen-output @@ -1,18 +1,18 @@ - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_bottom.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_west.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_top.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_east.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_east.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_south.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_west.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_bottom.0.txt. - Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_south.0.txt. + Loading Ascii data boundary file ASPECT_DIR/data/boundary-velocity/ascii-data/test/shell_3d_top.0.txt. Number of active cells: 24 (on 2 levels) Number of degrees of freedom: 1,277 (915+57+305) From e92cbedf2e05e222169a09e0f348143840cb2354 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 23 Oct 2024 11:26:34 +0200 Subject: [PATCH 4/7] Do not throw in destructor during unwinding --- include/aspect/plugins.h | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/include/aspect/plugins.h b/include/aspect/plugins.h index da4d97125f4..f1f5e40ed18 100644 --- a/include/aspect/plugins.h +++ b/include/aspect/plugins.h @@ -295,7 +295,16 @@ namespace aspect template ManagerBase::~ManagerBase() { - Assert (plugin_names.size() == plugin_objects.size(), ExcInternalError()); + // only check and throw if we are not unwinding the stack due + // to an active exception +#ifdef DEAL_II_HAVE_CXX17 + if (std::uncaught_exceptions() == 0) +#else + if (std::uncaught_exception() == false) +#endif + { + Assert (plugin_names.size() == plugin_objects.size(), ExcInternalError()); + } } From 56ef103e5e276b12fe843f61f07bd895ac2e8013 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 23 Oct 2024 11:27:06 +0200 Subject: [PATCH 5/7] Improve boundary check --- source/simulator/helper_functions.cc | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc index d54b5ee18a2..572a7b5c62e 100644 --- a/source/simulator/helper_functions.cc +++ b/source/simulator/helper_functions.cc @@ -2471,6 +2471,8 @@ namespace aspect "> is listed as having more " "than one type of temperature or heat flux boundary condition in the input file.")); + boundary_indicator_lists.emplace_back(boundary_composition_manager.get_fixed_composition_boundary_indicators()); + // Check that the periodic boundaries do not have other boundary conditions set using periodic_boundary_set = std::set, unsigned int>>; @@ -2516,8 +2518,11 @@ namespace aspect for (const auto &p : list) AssertThrow (all_boundary_indicators.find (p) != all_boundary_indicators.end(), - ExcMessage ("One of the boundary indicators listed in the input file " - "is not used by the geometry model.")); + ExcMessage ("Boundary indicator <" + + + Utilities::int_to_string(p) + + + "> is listed for a boundary condition, but is not used by the geometry model.")); if (parameters.nonlinear_solver == NonlinearSolver::single_Advection_no_Stokes) { @@ -2525,7 +2530,8 @@ namespace aspect for (unsigned int i=0; i<4; ++i) AssertThrow (boundary_indicator_lists[i].empty(), ExcMessage ("With the solver scheme `single Advection, no Stokes', " - "one cannot set boundary conditions for velocity.")); + "one cannot set boundary conditions for velocity or traction, " + "but a boundary condition has been set.")); } } From a9ddee7a3695536b908ba4532b777e88ee7f28ca Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 23 Oct 2024 13:14:39 +0200 Subject: [PATCH 6/7] Hide individual component masks --- include/aspect/boundary_traction/interface.h | 14 +++-- include/aspect/boundary_velocity/interface.h | 14 +++-- source/boundary_traction/interface.cc | 26 +++++++- source/boundary_velocity/interface.cc | 26 +++++++- source/simulator/assemblers/stokes.cc | 9 +-- source/simulator/assembly.cc | 4 +- source/simulator/core.cc | 28 +-------- source/simulator/helper_functions.cc | 66 +++++++++++--------- source/simulator/melt.cc | 7 +-- source/simulator/newton.cc | 2 +- source/simulator/stokes_matrix_free.cc | 13 ++-- 11 files changed, 117 insertions(+), 92 deletions(-) diff --git a/include/aspect/boundary_traction/interface.h b/include/aspect/boundary_traction/interface.h index a1235e0b073..ed9bcd8c460 100644 --- a/include/aspect/boundary_traction/interface.h +++ b/include/aspect/boundary_traction/interface.h @@ -136,13 +136,17 @@ namespace aspect get_active_plugin_boundary_indicators() const; /** - * Return a list of component masks that indicate for - * each active plugin which components it is responsible for. + * Return a component masks that indicates for + * each boundary which components are prescribed by + * this manager class. All plugins that are responsible + * for this boundary use the same component mask. * The list of plugin objects can be - * requested by calling get_active_plugins(). + * requested by calling get_active_plugins() and the + * list of boundaries they are responsible for is + * returnd by get_active_plugin_boundary_indicators(). */ - const std::vector & - get_active_plugin_component_masks() const; + ComponentMask + get_component_mask(const types::boundary_id boundary_id) const; /** * Declare the parameters of all known boundary traction plugins, as diff --git a/include/aspect/boundary_velocity/interface.h b/include/aspect/boundary_velocity/interface.h index 8f1ccf69fb5..9d38f415487 100644 --- a/include/aspect/boundary_velocity/interface.h +++ b/include/aspect/boundary_velocity/interface.h @@ -161,13 +161,17 @@ namespace aspect get_active_plugin_boundary_indicators() const; /** - * Return a list of component masks that indicate for - * each active plugin which components it is responsible for. + * Return a component masks that indicates for + * each boundary which components are prescribed by + * this manager class. All plugins that are responsible + * for this boundary use the same component mask. * The list of plugin objects can be - * requested by calling get_active_plugins(). + * requested by calling get_active_plugins() and the + * list of boundaries they are responsible for is + * returnd by get_active_plugin_boundary_indicators(). */ - const std::vector & - get_active_plugin_component_masks() const; + ComponentMask + get_component_mask(const types::boundary_id boundary_id) const; /** * Return a set of boundary indicators for which boundary diff --git a/source/boundary_traction/interface.cc b/source/boundary_traction/interface.cc index c234437feb5..32a99b2133f 100644 --- a/source/boundary_traction/interface.cc +++ b/source/boundary_traction/interface.cc @@ -132,10 +132,30 @@ namespace aspect template - const std::vector & - Manager::get_active_plugin_component_masks() const + ComponentMask + Manager::get_component_mask(const types::boundary_id boundary_id) const { - return component_masks; + Assert(prescribed_traction_boundary_indicators.find(boundary_id) != prescribed_traction_boundary_indicators.end(), + ExcMessage("The boundary traction manager class was asked for the " + "component mask of boundary indicator <" + + + Utilities::int_to_string(boundary_id) + + + "> with symbolic name <" + + + this->get_geometry_model().translate_id_to_symbol_name(boundary_id) + + + ">, but this boundary is not part of the active boundary traction plugins.")); + + // Since all component masks of plugins at the same boundary are identical, we can use + // the component mask of the first plugin we find that is responsible for this boundary. + for (unsigned int i=0; i - const std::vector & - Manager::get_active_plugin_component_masks() const + ComponentMask + Manager::get_component_mask(const types::boundary_id boundary_id) const { - return component_masks; + Assert(prescribed_velocity_boundary_indicators.find(boundary_id) != prescribed_velocity_boundary_indicators.end(), + ExcMessage("The boundary velocity manager class was asked for the " + "component mask of boundary indicator <" + + + Utilities::int_to_string(boundary_id) + + + "> with symbolic name <" + + + this->get_geometry_model().translate_id_to_symbol_name(boundary_id) + + + ">, but this boundary is not part of the active boundary velocity plugins.")); + + // Since all component masks of plugins at the same boundary are identical, we can use + // the component mask of the first plugin we find that is responsible for this boundary. + for (unsigned int i=0; i::face_iterator face = scratch.cell->face(scratch.face_number); - const auto &traction_bis = this->get_boundary_traction_manager().get_active_plugin_boundary_indicators(); + const auto &traction_bis = this->get_boundary_traction_manager().get_prescribed_boundary_traction_indicators(); - if (std::find(traction_bis.begin(), traction_bis.end(), face->boundary_id()) - != traction_bis.end()) + if (traction_bis.find(face->boundary_id()) != traction_bis.end()) { for (unsigned int q=0; qstokes_system_on_boundary_face.push_back( std::make_unique>()); @@ -767,7 +767,7 @@ namespace aspect = ( // see if we need to assemble traction boundary conditions. // only if so do we actually need to have an FEFaceValues object - !boundary_traction_manager.get_active_plugins().empty() + !boundary_traction_manager.get_prescribed_boundary_traction_indicators().empty() ? update_values | update_quadrature_points | diff --git a/source/simulator/core.cc b/source/simulator/core.cc index 515478a5f88..2c3eee98ce5 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -1376,32 +1376,6 @@ namespace aspect // for the prescribed velocity fields boundary_velocity_manager.update(); - // translate from component mask per plugin to component mask per boundary id - std::map boundary_component_masks; - for (unsigned int i=0; isecond == component_mask, - ExcMessage("Boundary indicator <" - + - Utilities::int_to_string(boundary_id) - + - "> with symbolic name <" - + - geometry_model->translate_id_to_symbol_name (boundary_id) - + - "> is listed as having different component masks " - "for different boundary velocity plugins. This is not allowed.")); - } - else - boundary_component_masks[boundary_id] = component_mask; - } - // put boundary conditions into constraints object for each boundary for (const auto boundary_id: boundary_velocity_manager.get_prescribed_boundary_velocity_indicators()) { @@ -1422,7 +1396,7 @@ namespace aspect boundary_id, vel, constraints, - boundary_component_masks.at(boundary_id)); + boundary_velocity_manager.get_component_mask(boundary_id)); } } diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc index 572a7b5c62e..4b30d67765e 100644 --- a/source/simulator/helper_functions.cc +++ b/source/simulator/helper_functions.cc @@ -2387,38 +2387,42 @@ namespace aspect boundary_indicator_lists.emplace_back(boundary_velocity_manager.get_prescribed_boundary_velocity_indicators()); boundary_indicator_lists.emplace_back(boundary_traction_manager.get_prescribed_boundary_traction_indicators()); - // Make sure that each combination of boundary velocity and boundary traction plugin + // Make sure that each combination of boundary velocity and boundary traction condition // either refers to different boundary indicators or to different components - const auto &velocity_boundary_ids = boundary_velocity_manager.get_active_plugin_boundary_indicators(); - const auto &traction_boundary_ids = boundary_traction_manager.get_active_plugin_boundary_indicators(); - - for (unsigned int i=0; i with symbolic name <" - + - geometry_model->translate_id_to_symbol_name (velocity_boundary_ids[i]) - + - "> is listed as having both " - "velocity and traction boundary conditions in the input file.")); - - // we have ensured the prescribed velocity and prescribed traction boundary conditions - // are compatible. In order to check them against the other boundary conditions, we - // need to remove the boundary indicator from one of the lists to make sure it only - // appears in one of them. We choose to remove the boundary indicator from the traction list. - boundary_indicator_lists[3].erase(traction_boundary_ids[j]); - } - } + for (const auto velocity_boundary_id: boundary_indicator_lists[2]) + { + bool found_compatible_duplicate_boundary_id = false; + for (const auto traction_boundary_id: boundary_indicator_lists[3]) + { + if (velocity_boundary_id == traction_boundary_id) + { + // if boundary ids are identical, make sure that the components are different + AssertThrow((boundary_velocity_manager.get_component_mask(velocity_boundary_id) & + boundary_traction_manager.get_component_mask(traction_boundary_id)) == + ComponentMask(introspection.n_components, false), + ExcMessage("Boundary indicator <" + + + Utilities::int_to_string(velocity_boundary_id) + + + "> with symbolic name <" + + + geometry_model->translate_id_to_symbol_name (velocity_boundary_id) + + + "> is listed as having both " + "velocity and traction boundary conditions in the input file.")); + + found_compatible_duplicate_boundary_id = true; + } + } + // we have ensured the prescribed velocity and prescribed traction boundary conditions + // for the current boundary id are compatible. In order to check them against the other + // boundary conditions, we need to remove the boundary indicator from one of the lists + // to make sure it only appears in one of them. We choose to remove the boundary + // indicator from the traction list. We cannot do that in the loop above, because it + // invalidates the range of the loop. + if (found_compatible_duplicate_boundary_id) + boundary_indicator_lists[3].erase(velocity_boundary_id); + } // for each combination of velocity boundary indicator lists, make sure that the // intersection is empty diff --git a/source/simulator/melt.cc b/source/simulator/melt.cc index e478f3d7a5c..4016c76b05e 100644 --- a/source/simulator/melt.cc +++ b/source/simulator/melt.cc @@ -1034,10 +1034,9 @@ namespace aspect Assert(face_no != numbers::invalid_unsigned_int,ExcInternalError()); const typename DoFHandler::face_iterator face = cell->face(face_no); - const auto &traction_bis = this->get_boundary_traction_manager().get_active_plugin_boundary_indicators(); + const auto &traction_bis = this->get_boundary_traction_manager().get_prescribed_boundary_traction_indicators(); - if (std::find(traction_bis.begin(), traction_bis.end(), face->boundary_id()) - != traction_bis.end()) + if (traction_bis.find(face->boundary_id()) != traction_bis.end()) { scratch.face_finite_element_values.reinit (cell, face_no); @@ -1689,7 +1688,7 @@ namespace aspect std::make_unique>()); // add the terms for traction boundary conditions - if (!this->get_boundary_traction_manager().get_active_plugins().empty()) + if (!this->get_boundary_traction_manager().get_prescribed_boundary_traction_indicators().empty()) { assemblers.stokes_system_on_boundary_face.push_back( std::make_unique> ()); diff --git a/source/simulator/newton.cc b/source/simulator/newton.cc index e22adb8d381..92efe142878 100644 --- a/source/simulator/newton.cc +++ b/source/simulator/newton.cc @@ -97,7 +97,7 @@ namespace aspect " defined that handles this formulation.")); // add the terms for traction boundary conditions - if (!this->get_boundary_traction_manager().get_active_plugins().empty()) + if (!this->get_boundary_traction_manager().get_prescribed_boundary_traction_indicators().empty()) { assemblers.stokes_system_on_boundary_face.push_back( std::make_unique>()); diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc index 5a716de77ce..825dff5c537 100644 --- a/source/simulator/stokes_matrix_free.cc +++ b/source/simulator/stokes_matrix_free.cc @@ -2487,24 +2487,23 @@ namespace aspect mg_constrained_dofs_A_block.initialize(dof_handler_v); std::set dirichlet_boundary = sim.boundary_velocity_manager.get_zero_boundary_velocity_indicators(); - for (unsigned int i=0; i Date: Wed, 23 Oct 2024 13:18:59 +0200 Subject: [PATCH 7/7] Fix typo --- include/aspect/boundary_traction/interface.h | 2 +- include/aspect/boundary_velocity/interface.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/include/aspect/boundary_traction/interface.h b/include/aspect/boundary_traction/interface.h index ed9bcd8c460..9cf00fb379e 100644 --- a/include/aspect/boundary_traction/interface.h +++ b/include/aspect/boundary_traction/interface.h @@ -143,7 +143,7 @@ namespace aspect * The list of plugin objects can be * requested by calling get_active_plugins() and the * list of boundaries they are responsible for is - * returnd by get_active_plugin_boundary_indicators(). + * returned by get_active_plugin_boundary_indicators(). */ ComponentMask get_component_mask(const types::boundary_id boundary_id) const; diff --git a/include/aspect/boundary_velocity/interface.h b/include/aspect/boundary_velocity/interface.h index 9d38f415487..88cf733f0e9 100644 --- a/include/aspect/boundary_velocity/interface.h +++ b/include/aspect/boundary_velocity/interface.h @@ -168,7 +168,7 @@ namespace aspect * The list of plugin objects can be * requested by calling get_active_plugins() and the * list of boundaries they are responsible for is - * returnd by get_active_plugin_boundary_indicators(). + * returned by get_active_plugin_boundary_indicators(). */ ComponentMask get_component_mask(const types::boundary_id boundary_id) const;