ASPECT
interface.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2022 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #ifndef _aspect_mesh_deformation_interface_h
23 #define _aspect_mesh_deformation_interface_h
24 
25 #include <aspect/plugins.h>
27 #include <aspect/global.h>
28 
29 #include <deal.II/fe/fe_system.h>
30 #include <deal.II/dofs/dof_handler.h>
31 #include <deal.II/lac/affine_constraints.h>
32 #include <deal.II/base/index_set.h>
33 #include <deal.II/base/mg_level_object.h>
34 #include <deal.II/lac/la_parallel_vector.h>
35 #include <deal.II/multigrid/mg_constrained_dofs.h>
36 #include <deal.II/multigrid/mg_transfer_matrix_free.h>
37 #include <deal.II/multigrid/mg_transfer_global_coarsening.templates.h>
39 
40 
41 namespace aspect
42 {
43  using namespace dealii;
44 
45  namespace Assemblers
46  {
53  template <int dim>
55  public SimulatorAccess<dim>
56  {
57  public:
58  ApplyStabilization(const double stabilization_theta);
59 
60  void
63 
64  private:
70  const double free_surface_theta;
71  };
72  }
73 
74  template <int dim> class Simulator;
75 
80  namespace MeshDeformation
81  {
89  template <int dim>
91  {
92  public:
96  virtual bool needs_surface_stabilization() const;
97 
98 
106  virtual
107  Tensor<1,dim>
108  compute_initial_deformation_on_boundary(const types::boundary_id boundary_indicator,
109  const Point<dim> &position) const;
110 
118  virtual
119  void
120  compute_velocity_constraints_on_boundary(const DoFHandler<dim> &mesh_deformation_dof_handler,
121  AffineConstraints<double> &mesh_velocity_constraints,
122  const std::set<types::boundary_id> &boundary_id) const;
123  };
124 
125 
126 
132  template <int dim>
133  class MeshDeformationHandler: public SimulatorAccess<dim>
134  {
135  public:
143 
147  ~MeshDeformationHandler() override;
148 
154  void initialize();
155 
160  void set_assemblers(const SimulatorAccess<dim> &simulator_access,
161  aspect::Assemblers::Manager<dim> &assemblers) const;
162 
167  void update();
168 
176  void execute();
177 
182  void setup_dofs();
183 
187  static
188  void declare_parameters (ParameterHandler &prm);
189 
193  void parse_parameters (ParameterHandler &prm);
194 
212  static
213  void
214  register_mesh_deformation
215  (const std::string &name,
216  const std::string &description,
217  void (*declare_parameters_function) (ParameterHandler &),
218  std::unique_ptr<Interface<dim>> (*factory_function) ());
219 
224  const std::map<types::boundary_id, std::vector<std::string>> &
225  get_active_mesh_deformation_names () const;
226 
231  const std::map<types::boundary_id,std::vector<std::unique_ptr<Interface<dim>>>> &
232  get_active_mesh_deformation_models () const;
233 
238  const std::set<types::boundary_id> &
239  get_active_mesh_deformation_boundary_indicators () const;
240 
245  const std::set<types::boundary_id> &
246  get_boundary_indicators_requiring_stabilization () const;
247 
253  const std::set<types::boundary_id> &
254  get_free_surface_boundary_indicators () const;
255 
259  double get_free_surface_theta () const;
260 
275  const LinearAlgebra::Vector &
276  get_initial_topography () const;
277 
282  const LinearAlgebra::Vector &
283  get_mesh_displacements () const;
284 
288  const DoFHandler<dim> &
289  get_mesh_deformation_dof_handler () const;
290 
300  template <typename MeshDeformationType,
301  typename = typename std::enable_if_t<std::is_base_of<Interface<dim>,MeshDeformationType>::value>>
302  bool
303  has_matching_mesh_deformation_object () const;
304 
316  template <typename MeshDeformationType,
317  typename = typename std::enable_if_t<std::is_base_of<Interface<dim>,MeshDeformationType>::value>>
318  const MeshDeformationType &
319  get_matching_mesh_deformation_object () const;
320 
321 
327  const Mapping<dim> &
328  get_level_mapping(const unsigned int level) const;
329 
339  static
340  void
341  write_plugin_graph (std::ostream &output_stream);
342 
346  DeclException1 (ExcMeshDeformationNameNotFound,
347  std::string,
348  << "Could not find entry <"
349  << arg1
350  << "> among the names of registered mesh deformation objects.");
351 
352  private:
360  void make_initial_constraints ();
361 
374  void make_constraints ();
375 
380  void compute_mesh_displacements ();
381 
386  void compute_mesh_displacements_gmg ();
387 
403  void set_initial_topography ();
404 
408  void interpolate_mesh_velocity ();
409 
413  void update_multilevel_deformation ();
414 
420 
425  const FESystem<dim> mesh_deformation_fe;
426 
431 
438 
446 
451 
461 
471 
476 
481 
486  AffineConstraints<double> mesh_velocity_constraints;
487 
492  AffineConstraints<double> mesh_vertex_constraints;
493 
498  std::map<types::boundary_id,std::vector<std::unique_ptr<Interface<dim>>>> mesh_deformation_objects;
499 
504  std::map<types::boundary_id, std::vector<std::string>> mesh_deformation_object_names;
505 
512 
518 
525  std::set<types::boundary_id> zero_mesh_deformation_boundary_indicators;
526 
532  std::set<types::boundary_id> free_surface_boundary_indicators;
533 
538  std::set<types::boundary_id> boundary_indicators_requiring_stabilization;
539 
541 
548 
552  MGLevelObject<std::unique_ptr<Mapping<dim>>> level_mappings;
553 
557  MGLevelObject<::LinearAlgebra::distributed::Vector<double>> level_displacements;
558 
563 
567  MGConstrainedDoFs mg_constrained_dofs;
568 
569  friend class Simulator<dim>;
570  friend class SimulatorAccess<dim>;
571  };
572 
573 
574 
575  template <int dim>
576  template <typename MeshDeformationType, typename>
577  inline
578  bool
580  {
581  for (const auto &object_iterator : mesh_deformation_objects)
582  for (const auto &p : object_iterator.second)
583  if (Plugins::plugin_type_matches<MeshDeformationType>(*p))
584  return true;
585 
586  return false;
587  }
588 
589 
590 
591  template <int dim>
592  template <typename MeshDeformationType, typename>
593  inline
594  const MeshDeformationType &
596  {
597  AssertThrow(has_matching_mesh_deformation_object<MeshDeformationType> (),
598  ExcMessage("You asked MeshDeformation::MeshDeformationHandler::get_matching_mesh_deformation_object() for a "
599  "mesh deformation object of type <" + boost::core::demangle(typeid(MeshDeformationType).name()) + "> "
600  "that could not be found in the current model. Activate this "
601  "mesh deformation in the input file."));
602 
603  for (const auto &object_iterator : mesh_deformation_objects)
604  for (const auto &p : object_iterator.second)
605  if (Plugins::plugin_type_matches<MeshDeformationType>(*p))
606  return Plugins::get_plugin_as_type<MeshDeformationType>(*p);
607 
608  // We will never get here, because we had the Assert above. Just to avoid warnings.
609  typename std::vector<std::unique_ptr<Interface<dim>>>::const_iterator mesh_def;
610  return Plugins::get_plugin_as_type<MeshDeformationType>(*(*mesh_def));
611 
612  }
613 
614 
621  template <int dim>
622  std::string
624 
625 
626 
634 #define ASPECT_REGISTER_MESH_DEFORMATION_MODEL(classname,name,description) \
635  template class classname<2>; \
636  template class classname<3>; \
637  namespace ASPECT_REGISTER_MESH_DEFORMATION_MODEL_ ## classname \
638  { \
639  aspect::internal::Plugins::RegisterHelper<aspect::MeshDeformation::Interface<2>,classname<2>> \
640  dummy_ ## classname ## _2d (&aspect::MeshDeformation::MeshDeformationHandler<2>::register_mesh_deformation, \
641  name, description); \
642  aspect::internal::Plugins::RegisterHelper<aspect::MeshDeformation::Interface<3>,classname<3>> \
643  dummy_ ## classname ## _3d (&aspect::MeshDeformation::MeshDeformationHandler<3>::register_mesh_deformation, \
644  name, description); \
645  }
646  }
647 }
648 
649 #endif
::TrilinosWrappers::MPI::BlockVector BlockVector
Definition: global.h:261
std::set< types::boundary_id > prescribed_mesh_deformation_boundary_indicators
Definition: interface.h:511
MGLevelObject<::LinearAlgebra::distributed::Vector< double > > level_displacements
Definition: interface.h:557
void write_plugin_graph(std::ostream &output_stream)
std::set< types::boundary_id > boundary_indicators_requiring_stabilization
Definition: interface.h:538
::TrilinosWrappers::MPI::Vector Vector
Definition: global.h:255
AffineConstraints< double > mesh_velocity_constraints
Definition: interface.h:486
std::set< types::boundary_id > tangential_mesh_deformation_boundary_indicators
Definition: interface.h:517
AffineConstraints< double > mesh_vertex_constraints
Definition: interface.h:492
MGLevelObject< std::unique_ptr< Mapping< dim > > > level_mappings
Definition: interface.h:552
std::map< types::boundary_id, std::vector< std::string > > mesh_deformation_object_names
Definition: interface.h:504
LinearAlgebra::BlockVector mesh_velocity
Definition: interface.h:437
std::map< types::boundary_id, std::vector< std::unique_ptr< Interface< dim > > > > mesh_deformation_objects
Definition: interface.h:498
MGTransferMF< dim, double > mg_transfer
Definition: interface.h:562
MGTransferMatrixFree< dim, NumberType > MGTransferMF
Definition: compat.h:45
std::string get_valid_model_names_pattern()
Definition: compat.h:59
void declare_parameters(ParameterHandler &prm)
Definition: compat.h:42
std::set< types::boundary_id > free_surface_boundary_indicators
Definition: interface.h:532
std::set< types::boundary_id > zero_mesh_deformation_boundary_indicators
Definition: interface.h:525
DeclException1(ProbabilityFunctionNegative, Point< dim >,<< "Your probability density function in the particle generator " "returned a negative probability density for the following position: "<< arg1<< ". Please check your function expression.")