ASPECT
interface.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2023 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 
234  // Forward declaration:
235  template <int dim>
237 
244  template <int dim>
246  {
247  public:
258  MaterialModelInputs(const unsigned int n_points,
259  const unsigned int n_comp);
260 
276  MaterialModelInputs(const DataPostprocessorInputs::Vector<dim> &input_data,
277  const Introspection<dim> &introspection,
278  const bool compute_strain_rate = true);
279 
280 
300  MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
301  const typename DoFHandler<dim>::active_cell_iterator &cell,
302  const Introspection<dim> &introspection,
303  const LinearAlgebra::BlockVector &solution_vector,
304  const bool compute_strain_rate = true);
305 
320 
324  MaterialModelInputs (MaterialModelInputs &&) noexcept = default;
325 
330  MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;
331 
335  MaterialModelInputs &operator= (MaterialModelInputs &&) = default;
336 
342  void reinit(const FEValuesBase<dim,dim> &fe_values,
343  const typename DoFHandler<dim>::active_cell_iterator &cell,
344  const Introspection<dim> &introspection,
345  const LinearAlgebra::BlockVector &solution_vector,
346  const bool compute_strain_rate = true);
347 
352  unsigned int n_evaluation_points() const;
353 
360  bool requests_property(const MaterialProperties::Property &property) const;
361 
366  std::vector<Point<dim>> position;
367 
371  std::vector<double> temperature;
372 
376  std::vector<double> pressure;
377 
382  std::vector<Tensor<1,dim>> pressure_gradient;
383 
391  std::vector<Tensor<1,dim>> velocity;
392 
398  std::vector<std::vector<double>> composition;
399 
412  std::vector<SymmetricTensor<2,dim>> strain_rate;
413 
429  typename DoFHandler<dim>::active_cell_iterator current_cell;
430 
439 
449  template <class AdditionalInputType>
450  AdditionalInputType *get_additional_input();
451 
456  template <class AdditionalInputType>
457  const AdditionalInputType *get_additional_input() const;
458 
464  std::vector<std::unique_ptr<AdditionalMaterialInputs<dim>>> additional_inputs;
465  };
466 
467 
468  // Forward declaration:
469  template <int dim>
471 
472 
479  template <int dim>
481  {
482  public:
493  MaterialModelOutputs (const unsigned int n_points,
494  const unsigned int n_comp);
495 
510 
514  MaterialModelOutputs (MaterialModelOutputs &&) noexcept = default;
515 
520  MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;
521 
525  MaterialModelOutputs &operator= (MaterialModelOutputs &&) noexcept = default;
526 
531  unsigned int n_evaluation_points() const;
532 
536  std::vector<double> viscosities;
537 
541  std::vector<double> densities;
542 
547  std::vector<double> thermal_expansion_coefficients;
548 
552  std::vector<double> specific_heat;
553 
557  std::vector<double> thermal_conductivities;
558 
563  std::vector<double> compressibilities;
564 
570  std::vector<double> entropy_derivative_pressure;
571 
578  std::vector<double> entropy_derivative_temperature;
579 
616  std::vector<std::vector<double>> reaction_terms;
617 
623  std::vector<std::unique_ptr<AdditionalMaterialOutputs<dim>>> additional_outputs;
624 
632  template <class AdditionalOutputType>
633  AdditionalOutputType *get_additional_output();
634 
638  template <class AdditionalOutputType>
639  const AdditionalOutputType *get_additional_output() const;
640 
645  void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
646  };
647 
648 
649 
664  namespace MaterialAveraging
665  {
705  {
717  };
718 
719 
726  std::string get_averaging_operation_names ();
727 
733  AveragingOperation parse_averaging_operation_name (const std::string &s);
734 
740  template <int dim>
741  void average (const AveragingOperation operation,
742  const typename DoFHandler<dim>::active_cell_iterator &cell,
743  const Quadrature<dim> &quadrature_formula,
744  const Mapping<dim> &mapping,
745  const MaterialProperties::Property &requested_properties,
746  MaterialModelOutputs<dim> &values_out);
747 
753  void average_property (const AveragingOperation operation,
754  const FullMatrix<double> &projection_matrix,
755  const FullMatrix<double> &expansion_matrix,
756  std::vector<double> &values_out);
757 
768  }
769 
788  template <int dim>
790  {
791  public:
796  virtual ~AdditionalMaterialInputs() = default;
797 
802  virtual void
803  fill (const LinearAlgebra::BlockVector &solution,
804  const FEValuesBase<dim> &fe_values,
805  const Introspection<dim> &introspection) = 0;
806  };
807 
808 
835  template <int dim>
837  {
838  public:
843  virtual ~AdditionalMaterialOutputs() = default;
844 
845  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
846  const FullMatrix<double> &/*projection_matrix*/,
847  const FullMatrix<double> &/*expansion_matrix*/)
848  {}
849  };
850 
851 
873  template <int dim>
875  {
876  public:
885  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
886 
897  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names,
898  const unsigned int n_points);
899 
903  ~NamedAdditionalMaterialOutputs() override;
904 
909  const std::vector<std::string> &get_names() const;
910 
915  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
916 
918  const FullMatrix<double> &/*projection_matrix*/,
919  const FullMatrix<double> &/*expansion_matrix*/) override
920  {}
921 
922 
927  std::vector<std::vector<double>> output_values;
928 
929  private:
930  const std::vector<std::string> names;
931  };
932 
933 
939  template <int dim>
941  {
942  public:
943  SeismicAdditionalOutputs(const unsigned int n_points);
944 
945  std::vector<double> get_nth_output(const unsigned int idx) const override;
946 
952  std::vector<double> vs;
953 
959  std::vector<double> vp;
960  };
961 
962 
983  template <int dim>
985  {
986  public:
987  ReactionRateOutputs (const unsigned int n_points,
988  const unsigned int n_comp);
989 
990  std::vector<double> get_nth_output(const unsigned int idx) const override;
991 
999  std::vector<std::vector<double>> reaction_rates;
1000  };
1001 
1002 
1003 
1009  template <int dim>
1011  {
1012  public:
1013  PhaseOutputs(const unsigned int n_points);
1014  };
1015 
1016 
1017 
1036  template <int dim>
1038  {
1039  public:
1040  PrescribedFieldOutputs (const unsigned int n_points,
1041  const unsigned int n_comp);
1042 
1043  std::vector<double> get_nth_output(const unsigned int idx) const override;
1044 
1052  std::vector<std::vector<double>> prescribed_field_outputs;
1053  };
1054 
1068  template <int dim>
1070  {
1071  public:
1072  PrescribedTemperatureOutputs (const unsigned int n_points);
1073 
1074  std::vector<double> get_nth_output(const unsigned int idx) const override;
1075 
1083  std::vector<double> prescribed_temperature_outputs;
1084  };
1085 
1092  template <int dim>
1094  {
1095  public:
1096  AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
1097  : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
1098  {}
1099 
1101  = default;
1102 
1104  const FullMatrix<double> &/*projection_matrix*/,
1105  const FullMatrix<double> &/*expansion_matrix*/) override
1106  {
1107  // TODO: not implemented
1108  }
1109 
1115  std::vector<Tensor<1,dim>> rhs_u;
1116 
1121  std::vector<double> rhs_p;
1122 
1127  std::vector<double> rhs_melt_pc;
1128  };
1129 
1130 
1131 
1147  template <int dim>
1149  {
1150  public:
1154  explicit PrescribedPlasticDilation (const unsigned int n_points);
1155 
1159  std::vector<double> get_nth_output(const unsigned int idx) const override;
1160 
1165  std::vector<double> dilation;
1166  };
1167 
1168 
1169 
1176  template <int dim>
1178  {
1179  public:
1180  ElasticOutputs(const unsigned int n_points)
1181  : elastic_force(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
1182  , viscoelastic_strain_rate(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
1183  {}
1184 
1185  ~ElasticOutputs() override
1186  = default;
1187 
1189  const FullMatrix<double> &/*projection_matrix*/,
1190  const FullMatrix<double> &/*expansion_matrix*/) override
1191  {
1192  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1193  return;
1194  }
1195 
1201  std::vector<SymmetricTensor<2,dim>> elastic_force;
1202 
1207  std::vector<SymmetricTensor<2,dim>> viscoelastic_strain_rate;
1208  };
1209 
1210 
1211 
1222  template <int dim>
1224  {
1225  public:
1226  EnthalpyOutputs(const unsigned int n_points)
1227  : enthalpies_of_fusion(n_points, numbers::signaling_nan<double>())
1228  {}
1229 
1230  virtual ~EnthalpyOutputs()
1231  = default;
1232 
1234  const FullMatrix<double> &/*projection_matrix*/,
1235  const FullMatrix<double> &/*expansion_matrix*/) override
1236  {
1237  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1238  return;
1239  }
1240 
1246  std::vector<double> enthalpies_of_fusion;
1247  };
1248 
1249 
1250 
1266  template <int dim>
1267  class Interface : public Plugins::InterfaceBase
1268  {
1269  public:
1277 
1285 
1298  get_model_dependence () const;
1299 
1309  virtual bool is_compressible () const = 0;
1318  virtual
1319  void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
1321 
1327  virtual
1328  void
1329  create_additional_named_outputs (MaterialModelOutputs &outputs) const;
1330 
1331 
1339  virtual
1340  void
1341  fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
1342  const LinearAlgebra::BlockVector &solution,
1343  const FEValuesBase<dim> &fe_values,
1344  const Introspection<dim> &introspection) const;
1345 
1346  protected:
1363  };
1364 
1380  template <int dim>
1381  void
1382  register_material_model (const std::string &name,
1383  const std::string &description,
1384  void (*declare_parameters_function) (ParameterHandler &),
1385  std::unique_ptr<Interface<dim>> (*factory_function) ());
1386 
1397  template <int dim>
1398  std::unique_ptr<Interface<dim>>
1399  create_material_model (const std::string &model_name);
1400 
1401 
1413  template <int dim>
1414  std::unique_ptr<Interface<dim>>
1415  create_material_model (ParameterHandler &prm);
1416 
1417 
1424  template <int dim>
1425  std::string
1427 
1428 
1434  template <int dim>
1435  void
1436  declare_parameters (ParameterHandler &prm);
1437 
1438 
1439 
1449  template <int dim>
1450  void
1451  write_plugin_graph (std::ostream &output_stream);
1452 
1453 
1454 
1455 // --------------------- template function definitions ----------------------------------
1456 
1457  template <int dim>
1458  template <class AdditionalInputType>
1460  {
1461  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1462  {
1463  AdditionalInputType *result = dynamic_cast<AdditionalInputType *> (additional_inputs[i].get());
1464  if (result)
1465  return result;
1466  }
1467  return nullptr;
1468  }
1469 
1470 
1471  template <int dim>
1472  template <class AdditionalInputType>
1473  const AdditionalInputType *MaterialModelInputs<dim>::get_additional_input() const
1474  {
1475  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1476  {
1477  const AdditionalInputType *result = dynamic_cast<const AdditionalInputType *> (additional_inputs[i].get());
1478  if (result)
1479  return result;
1480  }
1481  return nullptr;
1482  }
1483 
1484 
1485  template <int dim>
1486  template <class AdditionalOutputType>
1488  {
1489  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1490  {
1491  AdditionalOutputType *result = dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get());
1492  if (result)
1493  return result;
1494  }
1495  return nullptr;
1496  }
1497 
1498 
1499  template <int dim>
1500  template <class AdditionalOutputType>
1501  const AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output() const
1502  {
1503  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1504  {
1505  const AdditionalOutputType *result = dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get());
1506  if (result)
1507  return result;
1508  }
1509  return nullptr;
1510  }
1511 
1512 
1513  template <int dim>
1515  {
1516  Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
1517  this->additional_outputs = std::move(other.additional_outputs);
1518  }
1519 
1520 
1528 #define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
1529  template class classname<2>; \
1530  template class classname<3>; \
1531  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
1532  { \
1533  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2>> \
1534  dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
1535  name, description); \
1536  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3>> \
1537  dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
1538  name, description); \
1539  }
1540  }
1541 }
1542 
1543 
1544 #endif
void write_plugin_graph(std::ostream &output_stream)
AveragingOperation get_averaging_operation_for_viscosity(const AveragingOperation operation)
std::vector< double > entropy_derivative_pressure
Definition: interface.h:570
std::vector< double > entropy_derivative_temperature
Definition: interface.h:578
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:261
std::vector< std::vector< double > > composition
Definition: interface.h:398
std::vector< double > thermal_expansion_coefficients
Definition: interface.h:547
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:845
std::vector< Point< dim > > position
Definition: interface.h:366
std::vector< std::vector< double > > prescribed_field_outputs
Definition: interface.h:1052
AdditionalInputType * get_additional_input()
Definition: interface.h:1459
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< std::unique_ptr< AdditionalMaterialOutputs< dim > > > additional_outputs
Definition: interface.h:623
std::vector< double > thermal_conductivities
Definition: interface.h:557
ElasticOutputs(const unsigned int n_points)
Definition: interface.h:1180
std::vector< std::vector< double > > reaction_rates
Definition: interface.h:999
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:1103
std::vector< SymmetricTensor< 2, dim > > strain_rate
Definition: interface.h:412
std::vector< std::vector< double > > reaction_terms
Definition: interface.h:616
AdditionalOutputType * get_additional_output()
Definition: interface.h:1487
const std::vector< std::string > names
Definition: interface.h:930
Property operator|=(Property &d1, const Property d2)
Definition: interface.h:226
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1233
std::vector< Tensor< 1, dim > > velocity
Definition: interface.h:391
EnthalpyOutputs(const unsigned int n_points)
Definition: interface.h:1226
AveragingOperation parse_averaging_operation_name(const std::string &s)
Property operator|(const Property d1, const Property d2)
Definition: interface.h:220
AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
Definition: interface.h:1096
void move_additional_outputs_from(MaterialModelOutputs< dim > &other)
Definition: interface.h:1514
std::vector< double > enthalpies_of_fusion
Definition: interface.h:1246
std::vector< SymmetricTensor< 2, dim > > elastic_force
Definition: interface.h:1201
std::vector< std::vector< double > > output_values
Definition: interface.h:927
std::vector< std::unique_ptr< AdditionalMaterialInputs< dim > > > additional_inputs
Definition: interface.h:464
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1188
Definition: compat.h:59
DoFHandler< dim >::active_cell_iterator current_cell
Definition: interface.h:429
NonlinearDependence::ModelDependence model_dependence
Definition: interface.h:1362
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()
std::vector< Tensor< 1, dim > > pressure_gradient
Definition: interface.h:382
std::vector< SymmetricTensor< 2, dim > > viscoelastic_strain_rate
Definition: interface.h:1207
MaterialProperties::Property requested_properties
Definition: interface.h:438
Definition: compat.h:42
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:917