ASPECT
interface.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2024 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 #ifndef _aspect_material_model_interface_h
22 #define _aspect_material_model_interface_h
23 
24 #include <aspect/global.h>
25 #include <aspect/plugins.h>
27 
28 #include <deal.II/base/point.h>
29 #include <deal.II/base/quadrature.h>
30 #include <deal.II/base/symmetric_tensor.h>
31 #include <deal.II/base/parameter_handler.h>
32 #include <deal.II/dofs/dof_handler.h>
33 #include <deal.II/dofs/dof_accessor.h>
34 #include <deal.II/fe/mapping.h>
35 #include <deal.II/fe/fe_values.h>
36 #include <deal.II/fe/component_mask.h>
37 #include <deal.II/numerics/data_postprocessor.h>
38 #include <deal.II/base/signaling_nan.h>
39 
40 namespace aspect
41 {
42  template <int dim>
43  struct Introspection;
44 
45 
53  namespace MaterialModel
54  {
59  namespace NonlinearDependence
60  {
83  {
85 
86  none = 1,
88  pressure = 4,
91 
93  };
94 
95 
100  const Dependence d2)
101  {
102  return Dependence(static_cast<int>(d1) | static_cast<int>(d2));
103  }
104 
106  const Dependence d2)
107  {
108  d1 = (d1 | d2);
109  return d1;
110  }
111 
117  {
123 
129 
135 
141 
147 
153  ModelDependence ();
154  };
155 
164  bool
165  identifies_single_variable(const Dependence dependence);
166  }
167 
172  namespace MaterialProperties
173  {
183  enum Property
184  {
186 
187  none = 1,
189  density = 4,
199 
202  specific_heat |
207  viscosity |
212  };
213 
218  inline Property operator | (const Property d1,
219  const Property d2)
220  {
221  return Property(static_cast<int>(d1) | static_cast<int>(d2));
222  }
223 
225  const Property d2)
226  {
227  return (d1 | d2);
228  }
229  }
230 
231 
232  // Forward declaration:
233  template <int dim>
235 
242  template <int dim>
244  {
245  public:
256  MaterialModelInputs(const unsigned int n_points,
257  const unsigned int n_comp);
258 
274  MaterialModelInputs(const DataPostprocessorInputs::Vector<dim> &input_data,
275  const Introspection<dim> &introspection,
276  const bool compute_strain_rate = true);
277 
278 
298  MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
299  const typename DoFHandler<dim>::active_cell_iterator &cell,
300  const Introspection<dim> &introspection,
301  const LinearAlgebra::BlockVector &solution_vector,
302  const bool compute_strain_rate = true);
303 
318 
322  MaterialModelInputs (MaterialModelInputs &&) noexcept = default;
323 
328  MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;
329 
333  MaterialModelInputs &operator= (MaterialModelInputs &&) = default;
334 
340  void reinit(const FEValuesBase<dim,dim> &fe_values,
341  const typename DoFHandler<dim>::active_cell_iterator &cell,
342  const Introspection<dim> &introspection,
343  const LinearAlgebra::BlockVector &solution_vector,
344  const bool compute_strain_rate = true);
345 
350  unsigned int n_evaluation_points() const;
351 
358  bool requests_property(const MaterialProperties::Property &property) const;
359 
364  std::vector<Point<dim>> position;
365 
369  std::vector<double> temperature;
370 
374  std::vector<double> pressure;
375 
380  std::vector<Tensor<1,dim>> pressure_gradient;
381 
389  std::vector<Tensor<1,dim>> velocity;
390 
396  std::vector<std::vector<double>> composition;
397 
410  std::vector<SymmetricTensor<2,dim>> strain_rate;
411 
427  typename DoFHandler<dim>::active_cell_iterator current_cell;
428 
437 
446  template <class AdditionalInputType>
447  std::shared_ptr<AdditionalInputType>
448  get_additional_input_object();
449 
454  template <class AdditionalInputType>
455  std::shared_ptr<const AdditionalInputType>
456  get_additional_input_object() const;
457 
463  template <class AdditionalInputType>
464  DEAL_II_DEPRECATED
465  AdditionalInputType *
466  get_additional_input();
467 
473  template <class AdditionalInputType>
474  DEAL_II_DEPRECATED
475  const AdditionalInputType *
476  get_additional_input() const;
477 
487  template <class AdditionalInputType>
488  bool
489  has_additional_input_object() const;
490 
496  std::vector<std::shared_ptr<AdditionalMaterialInputs<dim>>> additional_inputs;
497  };
498 
499 
500  // Forward declaration:
501  template <int dim>
503 
504 
511  template <int dim>
513  {
514  public:
525  MaterialModelOutputs (const unsigned int n_points,
526  const unsigned int n_comp);
527 
542 
546  MaterialModelOutputs (MaterialModelOutputs &&) noexcept = default;
547 
552  MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;
553 
557  MaterialModelOutputs &operator= (MaterialModelOutputs &&) noexcept = default;
558 
563  unsigned int n_evaluation_points() const;
564 
568  std::vector<double> viscosities;
569 
573  std::vector<double> densities;
574 
579  std::vector<double> thermal_expansion_coefficients;
580 
584  std::vector<double> specific_heat;
585 
589  std::vector<double> thermal_conductivities;
590 
595  std::vector<double> compressibilities;
596 
602  std::vector<double> entropy_derivative_pressure;
603 
610  std::vector<double> entropy_derivative_temperature;
611 
648  std::vector<std::vector<double>> reaction_terms;
649 
655  std::vector<std::shared_ptr<AdditionalMaterialOutputs<dim>>> additional_outputs;
656 
666  template <class AdditionalOutputType>
667  std::shared_ptr<AdditionalOutputType>
668  get_additional_output_object();
669 
674  template <class AdditionalOutputType>
675  std::shared_ptr<const AdditionalOutputType>
676  get_additional_output_object() const;
677 
683  template <class AdditionalOutputType>
684  DEAL_II_DEPRECATED
685  AdditionalOutputType *
686  get_additional_output();
687 
693  template <class AdditionalOutputType>
694  DEAL_II_DEPRECATED
695  const AdditionalOutputType *
696  get_additional_output() const;
697 
707  template <class AdditionalInputType>
708  bool
709  has_additional_output_object() const;
710 
715  void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
716  };
717 
718 
719 
734  namespace MaterialAveraging
735  {
775  {
787  };
788 
789 
796  std::string get_averaging_operation_names ();
797 
803  AveragingOperation parse_averaging_operation_name (const std::string &s);
804 
810  template <int dim>
811  void average (const AveragingOperation operation,
812  const typename DoFHandler<dim>::active_cell_iterator &cell,
813  const Quadrature<dim> &quadrature_formula,
814  const Mapping<dim> &mapping,
815  const MaterialProperties::Property &requested_properties,
816  MaterialModelOutputs<dim> &values_out);
817 
823  void average_property (const AveragingOperation operation,
824  const FullMatrix<double> &projection_matrix,
825  const FullMatrix<double> &expansion_matrix,
826  std::vector<double> &values_out);
827 
838  }
839 
858  template <int dim>
860  {
861  public:
866  virtual ~AdditionalMaterialInputs() = default;
867 
872  virtual void
873  fill (const LinearAlgebra::BlockVector &solution,
874  const FEValuesBase<dim> &fe_values,
875  const Introspection<dim> &introspection) = 0;
876  };
877 
878 
905  template <int dim>
907  {
908  public:
913  virtual ~AdditionalMaterialOutputs() = default;
914 
915  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
916  const FullMatrix<double> &/*projection_matrix*/,
917  const FullMatrix<double> &/*expansion_matrix*/)
918  {}
919  };
920 
921 
943  template <int dim>
945  {
946  public:
955  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
956 
967  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names,
968  const unsigned int n_points);
969 
973  ~NamedAdditionalMaterialOutputs() override;
974 
979  const std::vector<std::string> &get_names() const;
980 
985  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
986 
988  const FullMatrix<double> &/*projection_matrix*/,
989  const FullMatrix<double> &/*expansion_matrix*/) override
990  {}
991 
992 
997  std::vector<std::vector<double>> output_values;
998 
999  private:
1000  const std::vector<std::string> names;
1001  };
1002 
1003 
1009  template <int dim>
1011  {
1012  public:
1013  SeismicAdditionalOutputs(const unsigned int n_points);
1014 
1015  std::vector<double> get_nth_output(const unsigned int idx) const override;
1016 
1022  std::vector<double> vs;
1023 
1029  std::vector<double> vp;
1030  };
1031 
1032 
1053  template <int dim>
1055  {
1056  public:
1057  ReactionRateOutputs (const unsigned int n_points,
1058  const unsigned int n_comp);
1059 
1060  std::vector<double> get_nth_output(const unsigned int idx) const override;
1061 
1069  std::vector<std::vector<double>> reaction_rates;
1070  };
1071 
1072 
1073 
1079  template <int dim>
1081  {
1082  public:
1083  PhaseOutputs(const unsigned int n_points);
1084  };
1085 
1086 
1087 
1106  template <int dim>
1108  {
1109  public:
1110  PrescribedFieldOutputs (const unsigned int n_points,
1111  const unsigned int n_comp);
1112 
1113  std::vector<double> get_nth_output(const unsigned int idx) const override;
1114 
1122  std::vector<std::vector<double>> prescribed_field_outputs;
1123  };
1124 
1138  template <int dim>
1140  {
1141  public:
1142  PrescribedTemperatureOutputs (const unsigned int n_points);
1143 
1144  std::vector<double> get_nth_output(const unsigned int idx) const override;
1145 
1153  std::vector<double> prescribed_temperature_outputs;
1154  };
1155 
1162  template <int dim>
1164  {
1165  public:
1166  AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
1167  : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
1168  {}
1169 
1171  = default;
1172 
1174  const FullMatrix<double> &/*projection_matrix*/,
1175  const FullMatrix<double> &/*expansion_matrix*/) override
1176  {
1177  // TODO: not implemented
1178  }
1179 
1185  std::vector<Tensor<1,dim>> rhs_u;
1186 
1191  std::vector<double> rhs_p;
1192 
1197  std::vector<double> rhs_melt_pc;
1198  };
1199 
1200 
1201 
1229  template <int dim>
1231  {
1232  public:
1236  explicit PrescribedPlasticDilation (const unsigned int n_points);
1237 
1241  std::vector<double> get_nth_output(const unsigned int idx) const override;
1242 
1247  std::vector<double> dilation_lhs_term;
1248 
1253  std::vector<double> dilation_rhs_term;
1254  };
1255 
1256 
1257 
1264  template <int dim>
1266  {
1267  public:
1268  ElasticOutputs(const unsigned int n_points)
1269  : elastic_force(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>()),
1270  viscoelastic_strain_rate(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
1271  {}
1272 
1273  ~ElasticOutputs() override
1274  = default;
1275 
1277  const FullMatrix<double> &/*projection_matrix*/,
1278  const FullMatrix<double> &/*expansion_matrix*/) override
1279  {
1280  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1281  return;
1282  }
1283 
1289  std::vector<SymmetricTensor<2,dim>> elastic_force;
1290 
1295  std::vector<SymmetricTensor<2,dim>> viscoelastic_strain_rate;
1296  };
1297 
1298 
1299 
1310  template <int dim>
1312  {
1313  public:
1314  EnthalpyOutputs(const unsigned int n_points)
1315  : enthalpies_of_fusion(n_points, numbers::signaling_nan<double>())
1316  {}
1317 
1318  virtual ~EnthalpyOutputs()
1319  = default;
1320 
1322  const FullMatrix<double> &/*projection_matrix*/,
1323  const FullMatrix<double> &/*expansion_matrix*/) override
1324  {
1325  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1326  return;
1327  }
1328 
1334  std::vector<double> enthalpies_of_fusion;
1335  };
1336 
1337 
1338 
1354  template <int dim>
1355  class Interface : public Plugins::InterfaceBase
1356  {
1357  public:
1365 
1373 
1386  get_model_dependence () const;
1387 
1397  virtual bool is_compressible () const = 0;
1406  virtual
1407  void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
1409 
1415  virtual
1416  void
1417  create_additional_named_outputs (MaterialModelOutputs &outputs) const;
1418 
1419 
1427  virtual
1428  void
1429  fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
1430  const LinearAlgebra::BlockVector &solution,
1431  const FEValuesBase<dim> &fe_values,
1432  const Introspection<dim> &introspection) const;
1433 
1434  protected:
1451  };
1452 
1468  template <int dim>
1469  void
1470  register_material_model (const std::string &name,
1471  const std::string &description,
1472  void (*declare_parameters_function) (ParameterHandler &),
1473  std::unique_ptr<Interface<dim>> (*factory_function) ());
1474 
1485  template <int dim>
1486  std::unique_ptr<Interface<dim>>
1487  create_material_model (const std::string &model_name);
1488 
1489 
1501  template <int dim>
1502  std::unique_ptr<Interface<dim>>
1503  create_material_model (ParameterHandler &prm);
1504 
1505 
1512  template <int dim>
1513  std::string
1515 
1516 
1522  template <int dim>
1523  void
1524  declare_parameters (ParameterHandler &prm);
1525 
1526 
1527 
1537  template <int dim>
1538  void
1539  write_plugin_graph (std::ostream &output_stream);
1540 
1541 
1542 
1543 // --------------------- template function definitions ----------------------------------
1544 
1545  template <int dim>
1546  template <class AdditionalInputType>
1547  std::shared_ptr<AdditionalInputType>
1549  {
1550  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1551  if (dynamic_cast<AdditionalInputType *> (additional_inputs[i].get()))
1552  return std::dynamic_pointer_cast<AdditionalInputType>(additional_inputs[i]);
1553 
1554  return nullptr;
1555  }
1556 
1557 
1558  template <int dim>
1559  template <class AdditionalInputType>
1560  std::shared_ptr<const AdditionalInputType>
1562  {
1563  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1564  if (dynamic_cast<AdditionalInputType *> (additional_inputs[i].get()))
1565  return std::dynamic_pointer_cast<const AdditionalInputType>(additional_inputs[i]);
1566 
1567  return nullptr;
1568  }
1569 
1570 
1571  template <int dim>
1572  template <class AdditionalOutputType>
1573  std::shared_ptr<AdditionalOutputType>
1575  {
1576  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1577  if (dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get()))
1578  return std::dynamic_pointer_cast<AdditionalOutputType>(additional_outputs[i]);
1579 
1580  return nullptr;
1581  }
1582 
1583 
1584  template <int dim>
1585  template <class AdditionalOutputType>
1586  std::shared_ptr<const AdditionalOutputType>
1588  {
1589  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1590  if (dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get()))
1591  return std::dynamic_pointer_cast<const AdditionalOutputType>(additional_outputs[i]);
1592 
1593  return nullptr;
1594  }
1595 
1596 
1597  // The following four functions are deprecated:
1598  template <int dim>
1599  template <class AdditionalInputType>
1600  AdditionalInputType *
1602  {
1603  return get_additional_input_object<AdditionalInputType>().get();
1604  }
1605 
1606 
1607  template <int dim>
1608  template <class AdditionalInputType>
1609  const AdditionalInputType *
1611  {
1612  return get_additional_input_object<AdditionalInputType>().get();
1613  }
1614 
1615 
1616  template <int dim>
1617  template <class AdditionalOutputType>
1618  AdditionalOutputType *
1620  {
1621  return get_additional_output_object<AdditionalOutputType>().get();
1622  }
1623 
1624 
1625  template <int dim>
1626  template <class AdditionalOutputType>
1627  const AdditionalOutputType *
1629  {
1630  return get_additional_output_object<AdditionalOutputType>().get();
1631  }
1632 
1633 
1634  template <int dim>
1635  template <class AdditionalInputType>
1636  bool
1638  {
1639  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1640  if (dynamic_cast<AdditionalInputType *> (additional_inputs[i].get()))
1641  return true;
1642 
1643  return false;
1644  }
1645 
1646 
1647  template <int dim>
1648  template <class AdditionalOutputType>
1649  bool
1651  {
1652  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1653  if (dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get()))
1654  return true;
1655 
1656  return false;
1657  }
1658 
1659 
1660 
1661  template <int dim>
1663  {
1664  Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
1665  this->additional_outputs = std::move(other.additional_outputs);
1666  }
1667 
1668 
1676 #define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
1677  template class classname<2>; \
1678  template class classname<3>; \
1679  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
1680  { \
1681  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2>> \
1682  dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
1683  name, description); \
1684  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3>> \
1685  dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
1686  name, description); \
1687  }
1688  }
1689 }
1690 
1691 
1692 #endif
void write_plugin_graph(std::ostream &output_stream)
AveragingOperation get_averaging_operation_for_viscosity(const AveragingOperation operation)
std::vector< std::shared_ptr< AdditionalMaterialInputs< dim > > > additional_inputs
Definition: interface.h:496
std::vector< double > entropy_derivative_pressure
Definition: interface.h:602
std::vector< double > entropy_derivative_temperature
Definition: interface.h:610
void average_property(const AveragingOperation operation, const FullMatrix< double > &projection_matrix, const FullMatrix< double > &expansion_matrix, std::vector< double > &values_out)
::TrilinosWrappers::MPI::BlockVector BlockVector
Definition: global.h:268
std::vector< std::vector< double > > composition
Definition: interface.h:396
std::vector< double > thermal_expansion_coefficients
Definition: interface.h:579
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:915
std::vector< Point< dim > > position
Definition: interface.h:364
std::vector< std::vector< double > > prescribed_field_outputs
Definition: interface.h:1122
void register_material_model(const std::string &name, const std::string &description, void(*declare_parameters_function)(ParameterHandler &), std::unique_ptr< Interface< dim >>(*factory_function)())
std::vector< double > thermal_conductivities
Definition: interface.h:589
ElasticOutputs(const unsigned int n_points)
Definition: interface.h:1268
std::vector< std::vector< double > > reaction_rates
Definition: interface.h:1069
void declare_parameters(ParameterHandler &prm)
bool identifies_single_variable(const Dependence dependence)
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1173
std::vector< SymmetricTensor< 2, dim > > strain_rate
Definition: interface.h:410
std::vector< std::vector< double > > reaction_terms
Definition: interface.h:648
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1321
std::shared_ptr< AdditionalOutputType > get_additional_output_object()
Definition: interface.h:1574
std::vector< Tensor< 1, dim > > velocity
Definition: interface.h:389
EnthalpyOutputs(const unsigned int n_points)
Definition: interface.h:1314
AveragingOperation parse_averaging_operation_name(const std::string &s)
std::shared_ptr< AdditionalInputType > get_additional_input_object()
Definition: interface.h:1548
AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
Definition: interface.h:1166
Dependence operator|=(Dependence &d1, const Dependence d2)
Definition: interface.h:105
void move_additional_outputs_from(MaterialModelOutputs< dim > &other)
Definition: interface.h:1662
std::vector< double > enthalpies_of_fusion
Definition: interface.h:1334
std::vector< SymmetricTensor< 2, dim > > elastic_force
Definition: interface.h:1289
DEAL_II_DEPRECATED AdditionalInputType * get_additional_input()
DEAL_II_DEPRECATED AdditionalOutputType * get_additional_output()
std::vector< std::vector< double > > output_values
Definition: interface.h:997
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1276
DoFHandler< dim >::active_cell_iterator current_cell
Definition: interface.h:427
NonlinearDependence::ModelDependence model_dependence
Definition: interface.h:1450
std::unique_ptr< Interface< dim > > create_material_model(const std::string &model_name)
void average(const AveragingOperation operation, const typename DoFHandler< dim >::active_cell_iterator &cell, const Quadrature< dim > &quadrature_formula, const Mapping< dim > &mapping, const MaterialProperties::Property &requested_properties, MaterialModelOutputs< dim > &values_out)
std::string get_valid_model_names_pattern()
Dependence operator|(const Dependence d1, const Dependence d2)
Definition: interface.h:99
std::vector< Tensor< 1, dim > > pressure_gradient
Definition: interface.h:380
std::vector< SymmetricTensor< 2, dim > > viscoelastic_strain_rate
Definition: interface.h:1295
MaterialProperties::Property requested_properties
Definition: interface.h:436
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:987
std::vector< std::shared_ptr< AdditionalMaterialOutputs< dim > > > additional_outputs
Definition: interface.h:655