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 #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  {
55  using namespace dealii;
56 
61  namespace NonlinearDependence
62  {
85  {
87 
88  none = 1,
90  pressure = 4,
93 
95  };
96 
97 
102  const Dependence d2)
103  {
104  return Dependence(static_cast<int>(d1) | static_cast<int>(d2));
105  }
106 
108  const Dependence d2)
109  {
110  d1 = (d1 | d2);
111  return d1;
112  }
113 
119  {
125 
131 
137 
143 
149 
155  ModelDependence ();
156  };
157 
166  bool
167  identifies_single_variable(const Dependence dependence);
168  }
169 
174  namespace MaterialProperties
175  {
185  enum Property
186  {
188 
189  none = 1,
191  density = 4,
201 
204  specific_heat |
209  viscosity |
214  };
215 
220  inline Property operator | (const Property d1,
221  const Property d2)
222  {
223  return Property(static_cast<int>(d1) | static_cast<int>(d2));
224  }
225 
227  const Property d2)
228  {
229  return (d1 | d2);
230  }
231  }
232 
233  template <int dim> class AdditionalMaterialInputs;
234 
241  template <int dim>
243  {
254  MaterialModelInputs(const unsigned int n_points,
255  const unsigned int n_comp);
256 
272  MaterialModelInputs(const DataPostprocessorInputs::Vector<dim> &input_data,
273  const Introspection<dim> &introspection,
274  const bool compute_strain_rate = true);
275 
276 
296  MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
297  const typename DoFHandler<dim>::active_cell_iterator &cell,
298  const Introspection<dim> &introspection,
299  const LinearAlgebra::BlockVector &solution_vector,
300  const bool compute_strain_rate = true);
301 
316 
320  MaterialModelInputs (MaterialModelInputs &&) noexcept = default;
321 
326  MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;
327 
331  MaterialModelInputs &operator= (MaterialModelInputs &&) = default;
332 
338  void reinit(const FEValuesBase<dim,dim> &fe_values,
339  const typename DoFHandler<dim>::active_cell_iterator &cell,
340  const Introspection<dim> &introspection,
341  const LinearAlgebra::BlockVector &solution_vector,
342  const bool compute_strain_rate = true);
343 
348  unsigned int n_evaluation_points() const;
349 
356  bool requests_property(const MaterialProperties::Property &property) const;
357 
362  std::vector<Point<dim>> position;
363 
367  std::vector<double> temperature;
368 
372  std::vector<double> pressure;
373 
378  std::vector<Tensor<1,dim>> pressure_gradient;
379 
387  std::vector<Tensor<1,dim>> velocity;
388 
394  std::vector<std::vector<double>> composition;
395 
408  std::vector<SymmetricTensor<2,dim>> strain_rate;
409 
425  typename DoFHandler<dim>::active_cell_iterator current_cell;
426 
435 
441  std::vector<std::unique_ptr<AdditionalMaterialInputs<dim>>> additional_inputs;
442 
449  template <class AdditionalInputType>
450  AdditionalInputType *get_additional_input();
451 
455  template <class AdditionalInputType>
456  const AdditionalInputType *get_additional_input() const;
457  };
458 
459 
460  template <int dim> class AdditionalMaterialOutputs;
461 
462 
469  template <int dim>
471  {
482  MaterialModelOutputs (const unsigned int n_points,
483  const unsigned int n_comp);
484 
499 
503  MaterialModelOutputs (MaterialModelOutputs &&) noexcept = default;
504 
509  MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;
510 
514  MaterialModelOutputs &operator= (MaterialModelOutputs &&) noexcept = default;
515 
520  unsigned int n_evaluation_points() const;
521 
525  std::vector<double> viscosities;
526 
530  std::vector<double> densities;
531 
536  std::vector<double> thermal_expansion_coefficients;
537 
541  std::vector<double> specific_heat;
542 
546  std::vector<double> thermal_conductivities;
547 
552  std::vector<double> compressibilities;
553 
559  std::vector<double> entropy_derivative_pressure;
560 
567  std::vector<double> entropy_derivative_temperature;
568 
605  std::vector<std::vector<double>> reaction_terms;
606 
612  std::vector<std::unique_ptr<AdditionalMaterialOutputs<dim>>> additional_outputs;
613 
621  template <class AdditionalOutputType>
622  AdditionalOutputType *get_additional_output();
623 
627  template <class AdditionalOutputType>
628  const AdditionalOutputType *get_additional_output() const;
629 
634  void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
635  };
636 
637 
638 
653  namespace MaterialAveraging
654  {
694  {
705  };
706 
707 
714  std::string get_averaging_operation_names ();
715 
721  AveragingOperation parse_averaging_operation_name (const std::string &s);
722 
728  template <int dim>
729  void average (const AveragingOperation operation,
730  const typename DoFHandler<dim>::active_cell_iterator &cell,
731  const Quadrature<dim> &quadrature_formula,
732  const Mapping<dim> &mapping,
733  MaterialModelOutputs<dim> &values_out);
734 
740  void average_property (const AveragingOperation operation,
741  const FullMatrix<double> &projection_matrix,
742  const FullMatrix<double> &expansion_matrix,
743  std::vector<double> &values_out);
744  }
745 
746 
765  template <int dim>
767  {
768  public:
769  virtual ~AdditionalMaterialInputs() = default;
770 
775  virtual void
776  fill (const LinearAlgebra::BlockVector &solution,
777  const FEValuesBase<dim> &fe_values,
778  const Introspection<dim> &introspection) = 0;
779  };
780 
781 
808  template <int dim>
810  {
811  public:
816  virtual ~AdditionalMaterialOutputs() = default;
817 
818  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
819  const FullMatrix<double> &/*projection_matrix*/,
820  const FullMatrix<double> &/*expansion_matrix*/)
821  {}
822  };
823 
824 
846  template <int dim>
848  {
849  public:
858  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
859 
870  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names,
871  const unsigned int n_points);
872 
876  ~NamedAdditionalMaterialOutputs() override;
877 
882  const std::vector<std::string> &get_names() const;
883 
888  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
889 
891  const FullMatrix<double> &/*projection_matrix*/,
892  const FullMatrix<double> &/*expansion_matrix*/) override
893  {}
894 
895 
900  std::vector<std::vector<double>> output_values;
901 
902  private:
903  const std::vector<std::string> names;
904  };
905 
906 
912  template <int dim>
914  {
915  public:
916  SeismicAdditionalOutputs(const unsigned int n_points);
917 
918  std::vector<double> get_nth_output(const unsigned int idx) const override;
919 
925  std::vector<double> vs;
926 
932  std::vector<double> vp;
933  };
934 
935 
956  template <int dim>
958  {
959  public:
960  ReactionRateOutputs (const unsigned int n_points,
961  const unsigned int n_comp);
962 
963  std::vector<double> get_nth_output(const unsigned int idx) const override;
964 
972  std::vector<std::vector<double>> reaction_rates;
973  };
974 
975 
976 
982  template <int dim>
984  {
985  public:
986  PhaseOutputs(const unsigned int n_points);
987  };
988 
989 
990 
1009  template <int dim>
1011  {
1012  public:
1013  PrescribedFieldOutputs (const unsigned int n_points,
1014  const unsigned int n_comp);
1015 
1016  std::vector<double> get_nth_output(const unsigned int idx) const override;
1017 
1025  std::vector<std::vector<double>> prescribed_field_outputs;
1026  };
1027 
1041  template <int dim>
1043  {
1044  public:
1045  PrescribedTemperatureOutputs (const unsigned int n_points);
1046 
1047  std::vector<double> get_nth_output(const unsigned int idx) const override;
1048 
1056  std::vector<double> prescribed_temperature_outputs;
1057  };
1058 
1065  template <int dim>
1067  {
1068  public:
1069  AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
1070  : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
1071  {}
1072 
1074  = default;
1075 
1077  const FullMatrix<double> &/*projection_matrix*/,
1078  const FullMatrix<double> &/*expansion_matrix*/) override
1079  {
1080  // TODO: not implemented
1081  }
1082 
1088  std::vector<Tensor<1,dim>> rhs_u;
1089 
1094  std::vector<double> rhs_p;
1095 
1100  std::vector<double> rhs_melt_pc;
1101  };
1102 
1103 
1104 
1120  template <int dim>
1122  {
1123  public:
1127  explicit PrescribedPlasticDilation (const unsigned int n_points);
1128 
1132  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
1133 
1138  std::vector<double> dilation;
1139  };
1140 
1141 
1142 
1149  template <int dim>
1151  {
1152  public:
1153  ElasticOutputs(const unsigned int n_points)
1154  : elastic_force(n_points, numbers::signaling_nan<Tensor<2,dim>>() )
1155  {}
1156 
1157  ~ElasticOutputs() override
1158  = default;
1159 
1161  const FullMatrix<double> &/*projection_matrix*/,
1162  const FullMatrix<double> &/*expansion_matrix*/) override
1163  {
1164  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1165  return;
1166  }
1167 
1173  std::vector<Tensor<2,dim>> elastic_force;
1174  };
1175 
1176 
1177 
1188  template <int dim>
1190  {
1191  public:
1192  EnthalpyOutputs(const unsigned int n_points)
1193  : enthalpies_of_fusion(n_points, numbers::signaling_nan<double>())
1194  {}
1195 
1196  virtual ~EnthalpyOutputs()
1197  = default;
1198 
1199  virtual void average (const MaterialAveraging::AveragingOperation operation,
1200  const FullMatrix<double> &/*projection_matrix*/,
1201  const FullMatrix<double> &/*expansion_matrix*/)
1202  {
1203  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1204  return;
1205  }
1206 
1212  std::vector<double> enthalpies_of_fusion;
1213  };
1214 
1215 
1216 
1232  template <int dim>
1233  class Interface
1234  {
1235  public:
1250 
1255  virtual ~Interface() = default;
1256 
1262  virtual
1263  void
1264  initialize ();
1265 
1270  virtual void update ();
1271 
1284  get_model_dependence () const;
1285 
1295  virtual bool is_compressible () const = 0;
1304  virtual
1305  void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
1317  static
1318  void
1319  declare_parameters (ParameterHandler &prm);
1320 
1327  virtual
1328  void
1329  parse_parameters (ParameterHandler &prm);
1339  virtual
1340  void
1341  create_additional_named_outputs (MaterialModelOutputs &outputs) const;
1342 
1343 
1351  virtual
1352  void
1353  fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
1354  const LinearAlgebra::BlockVector &solution,
1355  const FEValuesBase<dim> &fe_values,
1356  const Introspection<dim> &introspection) const;
1357 
1358  protected:
1375  };
1376 
1392  template <int dim>
1393  void
1394  register_material_model (const std::string &name,
1395  const std::string &description,
1396  void (*declare_parameters_function) (ParameterHandler &),
1397  std::unique_ptr<Interface<dim>> (*factory_function) ());
1398 
1409  template <int dim>
1410  std::unique_ptr<Interface<dim>>
1411  create_material_model (const std::string &model_name);
1412 
1413 
1425  template <int dim>
1426  std::unique_ptr<Interface<dim>>
1427  create_material_model (ParameterHandler &prm);
1428 
1429 
1436  template <int dim>
1437  std::string
1439 
1440 
1446  template <int dim>
1447  void
1448  declare_parameters (ParameterHandler &prm);
1449 
1450 
1451 
1461  template <int dim>
1462  void
1463  write_plugin_graph (std::ostream &output_stream);
1464 
1465 
1466 
1467 // --------------------- template function definitions ----------------------------------
1468 
1469  template <int dim>
1470  template <class AdditionalInputType>
1472  {
1473  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1474  {
1475  AdditionalInputType *result = dynamic_cast<AdditionalInputType *> (additional_inputs[i].get());
1476  if (result)
1477  return result;
1478  }
1479  return nullptr;
1480  }
1481 
1482 
1483  template <int dim>
1484  template <class AdditionalInputType>
1485  const AdditionalInputType *MaterialModelInputs<dim>::get_additional_input() const
1486  {
1487  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1488  {
1489  const AdditionalInputType *result = dynamic_cast<const AdditionalInputType *> (additional_inputs[i].get());
1490  if (result)
1491  return result;
1492  }
1493  return nullptr;
1494  }
1495 
1496 
1497  template <int dim>
1498  template <class AdditionalOutputType>
1500  {
1501  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1502  {
1503  AdditionalOutputType *result = dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get());
1504  if (result)
1505  return result;
1506  }
1507  return nullptr;
1508  }
1509 
1510 
1511  template <int dim>
1512  template <class AdditionalOutputType>
1513  const AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output() const
1514  {
1515  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1516  {
1517  const AdditionalOutputType *result = dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get());
1518  if (result)
1519  return result;
1520  }
1521  return nullptr;
1522  }
1523 
1524 
1525  template <int dim>
1527  {
1528  Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
1529  this->additional_outputs = std::move(other.additional_outputs);
1530  }
1531 
1532 
1540 #define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
1541  template class classname<2>; \
1542  template class classname<3>; \
1543  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
1544  { \
1545  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2>> \
1546  dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
1547  name, description); \
1548  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3>> \
1549  dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
1550  name, description); \
1551  }
1552  }
1553 }
1554 
1555 
1556 #endif
void write_plugin_graph(std::ostream &output_stream)
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:229
AdditionalInputType * get_additional_input()
Definition: interface.h:1471
std::vector< SymmetricTensor< 2, dim > > strain_rate
Definition: interface.h:408
std::vector< double > thermal_expansion_coefficients
Definition: interface.h:536
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:818
MaterialProperties::Property requested_properties
Definition: interface.h:434
std::vector< std::vector< double > > prescribed_field_outputs
Definition: interface.h:1025
std::vector< std::vector< double > > composition
Definition: interface.h:394
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< Point< dim > > position
Definition: interface.h:362
ElasticOutputs(const unsigned int n_points)
Definition: interface.h:1153
AdditionalOutputType * get_additional_output()
Definition: interface.h:1499
std::vector< std::vector< double > > reaction_rates
Definition: interface.h:972
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:1076
void average(const AveragingOperation operation, const typename DoFHandler< dim >::active_cell_iterator &cell, const Quadrature< dim > &quadrature_formula, const Mapping< dim > &mapping, MaterialModelOutputs< dim > &values_out)
const std::vector< std::string > names
Definition: interface.h:903
Property operator|=(Property &d1, const Property d2)
Definition: interface.h:226
std::vector< std::vector< double > > reaction_terms
Definition: interface.h:605
EnthalpyOutputs(const unsigned int n_points)
Definition: interface.h:1192
AveragingOperation parse_averaging_operation_name(const std::string &s)
DoFHandler< dim >::active_cell_iterator current_cell
Definition: interface.h:425
Property operator|(const Property d1, const Property d2)
Definition: interface.h:220
AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
Definition: interface.h:1069
std::vector< double > enthalpies_of_fusion
Definition: interface.h:1212
std::vector< std::unique_ptr< AdditionalMaterialInputs< dim > > > additional_inputs
Definition: interface.h:441
void move_additional_outputs_from(MaterialModelOutputs< dim > &other)
Definition: interface.h:1526
std::vector< Tensor< 1, dim > > pressure_gradient
Definition: interface.h:378
std::vector< std::vector< double > > output_values
Definition: interface.h:900
std::vector< Tensor< 2, dim > > elastic_force
Definition: interface.h:1173
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1160
std::vector< double > entropy_derivative_pressure
Definition: interface.h:559
std::vector< double > thermal_conductivities
Definition: interface.h:546
std::vector< double > entropy_derivative_temperature
Definition: interface.h:567
NonlinearDependence::ModelDependence model_dependence
Definition: interface.h:1374
std::unique_ptr< Interface< dim > > create_material_model(const std::string &model_name)
std::string get_valid_model_names_pattern()
std::vector< Tensor< 1, dim > > velocity
Definition: interface.h:387
virtual void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:1199
Definition: compat.h:88
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:890
std::vector< std::unique_ptr< AdditionalMaterialOutputs< dim > > > additional_outputs
Definition: interface.h:612